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7iBSTRACT 


The  second-order  solution  of  the  problem  of  the  inter- 
action of  a  train  or  regular  waves  with  a  completely  sub- 
merged, horizontal  circular  cylinder  in  finite  depth  water 
is  presented  for  the  two-dimensional  case.   The  incident 
wave  is  developed  as  a  second-order  Stokes1  wave  by  use  of 
a  perturbation  method  and  the  solution  of  both  the  first- 
order  and  second-order  scattering  potentials  is  obtained 
numerically  using  the  Green's  function  approach.   The  hydro- 
dynamic  pressure  and  resulting  first-order  and  second-order 
force  coefficients  are  determined  numerically  and  presented 
for  various  values  of  water  depth,  cylinder  depth  of 
submergence,  and  wave  length. 
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I.   INTRODUCTION 

The  solution  to  the  two-dimensional  linear  wave/structure 
interaction  problem  has  been  well-studied  for  both  the  finite 
and  infinite  depth  cases.   The  fluid  motion  resulting  from 
the  interaction  of  a  train  of  regular  waves  with  a  submerged 
horizontal  circular  cylinder  in  infinite  depth  water  was 
first  studied  by  Dean  [Ref.  1] .   Ursell  [Ref.  7]  later  studied 
the  problem  anew  placing  the  solution  on  a  rigorous  basis. 
More  recently,  Ogilvie  [Ref.  6]  reconsidered  the  problem.   He 
computed  the  first-order  oscillatory  forces  and  second-order 
steady-state  forces  for  the  following  cases:   (a)  the  cylinder 
held  fixed,  (b)  the  cylinder  forced  to  oscillate  sinusoidally 
in  still  water,  and  (c)  the  neutrally  buoyant  cylinder  allowed 
to  respond  to  wave  motion. 

It  can  easily  be  shown  that  the  steady-state  part  of  the 
second-order  forces  can  be  computed  from  the  first-order 
potential  alone.   Accordingly,  Ogilvie  did  not  obtain  a 
complete  second-order  solution;  the  steady-state  forces 
arise  from  the  first-order  velocity  squared  terms  in  Ber- 
noulli's equation.   In  addition  to  the  steady-state  forces, 
second-order  oscillatory  forces  also  exist  which  have  a 
frequency  twice  that  of  the  first-order  forces  and  represent 
the  subject  of  the  present  work. 

This  thesis  reconsiders  again  the  submerged  horizontal 
cylinder  problem  but  the  extension  is  made  to  include  water 
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of  finite  depth.   The  cylinder  is  considered  to  be  completely 
submerged  and  held  fixed  in  a  train  of  regular  waves.   The 
objective  is  to  determine  the  complete  second-order  solution 
and,  accordingly,  determine  the  oscillatory  second-order 
forces  as  well  as  the  steady-state  second-order  forces. 

The  problem  is  treated  as  a  regular  perturbation  problem 
in  the  small  parameter,  incident  wave  height/cylinder  radius. 
Proceeding  in  this  way  the  incident  wave  appears  as  a  second- 
order  Stoke' s  wave.   The  boundary-value  problems  for  both 
the  first-order  and  second-order  potentials  are  established 
and  the  solutions  to  both  are  obtained  by  use  of  the  Green's 
function  method. 
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II.   THEORETICAL  DEVELOPMENT 

A.   FORMULATION  OF  THE  PROBLEM 

A  systematic  representation  of  the  problem  considered  is 
illustrated  in  Fig.  (1) .   A  rigid  right  cylinder  of  radius  a, 
submerged  to  a  depth  d  in  water  of  depth  h,  is  subjected  to 
a  train  of  regular  waves  propagating  in  the  positive  x  direc- 
tion.  The  basic  problem  is  that  of  calculating  the  hydro- 
dynamic  pressure  on  the  immersed  cylinder  and  the  resulting 
force  correct  to  the  second-order  in  wave  height  of  the 
incident  wave. 

Assuming  the  fluid  to  be  irrotational ,  a  velocity 
potential,  $,  may  be  defined  as: 

q  =  V$(x,y,t)  (1) 

where  q  denotes  the  fluid  velocity  vector.   (The  barred 
quantities  denote  dimensional  quantities.)   Moreover,  assuming 
the  fluid  to  be  incompressible,  the  velocity  potential  must 
satisfy  the  Laplace  equation, 


V2$(x,y,t)  =  0  (2) 


within  the  fluid  region. 

In  addition  to  the  differential  equation,  Eq.  (2) ,  $ 
must  satisfy  certain  boundary  conditions.   In  specific, 
these  are: 
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1.  The  z.ero  normal  velocity  condition  on  the  bottom, 
defined  by  y  =  -  h,  where  H"  denotes  the  mean  fluid 
depth. 

2.  The  zero  normal  velocity  condition  on  the  surface, 
Sn(x,y)  =  0,  of  the  cylinder. 

3.  The  kinematic  boundary  condition  on  the  free  surface, 
S2 (x,y,t)  =  n(x,t)  =  0. 

4.  The  dynamic  boundary  condition  on  the  free  surface. 
These  boundary  conditions  may  be  stated  mathematically  as 
follows: 


$-(x,-h,t)  =  0  (3) 


V$-VS1  =  0  (4) 


dS2 
q*vs2  +  JT-=   °  (5) 


«^(x,ri/t)  +  I[<j>_(x,n,t)2  +  0-(x,n,t)2]  +  gn(x,t)  =  gK 

(6) 

where  n  denotes  the  elevation  of  the  free  surface,  g  denotes 
the  acceleration  of  gravity,  and  K  denotes  the  Bernoulli 
constant. 

For  convenience  in  carrying  out  the  solution,  the 
variables  are  next  made  dimensionless  using  the  cylinder 
radius,  and  the  wave  frequency,  o,  as  follows: 
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x  =  x/a      y  =  y/a      d  =  d/a      h  =  h/a 

H  =  H/2a     K  =  K/a      v  =  a2a/g    t  =  at        (7) 


<J>  =  a$/ga    ri  =  n/a 

H  denotes  the  dimensionless  wave  height,  K  denotes  the 
dimensionless  Bernoulli  constant,  and  t  denotes  the 
dimensionless  time. 

Utilizing  the  parameters  defined  in  Eq.  (7),  Eqs .  (2-6) 
may  be  rewritten  to  concisely  define  the  boundary-value 
problem  in  dimensionless  form  as: 


2 

V    (j>(x,y,t)    =    0  in    the    fluid   region  (8) 


A    (x,-h,t)    =   0  (9) 


<j>n(x,y,t)    =0  on   S-^Xjy)    =   0  (10) 


4>x(x,n,t)nx(x,t)    -  <j>   (x,n,t)   +  vnt(x,t)   =  o  (11) 


At(x,n,t)   +  j[$x(xtr\,t)2  +  a   (x,n,t)2]   +  n(x,t)   =  k     (12) 


B.   PERTURBATION  EXPANSION 

Expanding  <j>  and  n  in  terms  of  a  perturbation  parameter, 
e,  provides  a  means  for  converting  the  nonlinear  boundary-value 
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problem  into  a  series  of  linear  problems.  The  solution  to 
each  linear  problem  may  be  tractable  whereas  the  nonlinear 
problem  is  unsolvable. 

Since  <f>,  n,  and  K  are  functions  of  the  small  parameter, 
e,  they  may  be  written  in  the  power  series: 


and 


<|>(x,y,t;e)  =  £  en<J>  (x,y ,  t)  (13) 

n=l    n 


n(x/t;e)  =  £  enn  (x,t)  (14) 

n=l    n 


K  -  t2K2  +  0(£3)  (15) 


In  Eqs.  (11)  and  (12),  $  contains  n,  and,  therefore,  e 
implicitly.   This  can  be  converted  to  an  explicit  form  in 
e  by  use  of  the  Taylor  series  expansion: 


* 


oo  mm 

(x,n,t)  =  £  "(x't;c)  '  [5  4>(*'y't>]  _n  (i6) 


m=0    ml       3y 


m       Jy=0 


The  perturbation  parameter,  e,  is  related  to  the  wave  height 
in  an,  as  yet,  unknown  manner. 

Equations  (13-16)  as  well  as  appropriate  derivatives 
similar  to  Eqs.  (13-16)  may  now  be  substituted  into  the 
boundary-value  problem  given  in  Eqs.  (8-12)  to  obtain  a 
series  of  linear  boundary-value  problems  for  $. ,  $-,  etc. 
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That  is,  upon  equating  coefficients  of  like  powers  of  z, 
the  following  problems  for  the  first  two  terms  in  the 
expansion  for  <J>  can  be  precipitated: 
First-order  (e) : 


V2(J)1(x/y,t)    =   0  (17) 


<f>ly(x,~h,t)    =    0  (18) 

4>ln(x,y,t)    =0               on  S1(x/y)  =   0                            (19) 

<j>ly(x,0,t)    -   vnlt(x/t)  =   0  (20) 

n1(x/t)   +  <j>lt(x,o,t)   =  0  (21) 


2 

Second-order    (e    ) : 


V2<j>2(x;y,t)    =   0  (22) 


<J>2    (x,-h,t)    =   0  (23) 

<f>2n(x,y,t)    =0  on   S]L(x,y)    =   0  (24) 


<J>2 y(x,0,t)    -  vn2t(x/t) 


-n1(x,t)<J>1      (x,0,t)    +   cf>lx(x,o,t)nlx(x,t) 


(25) 


n2(x,t)  +  <j>2t(x,o,t)  -  -n1(x,t)<^1  t(x,o,t)   - 
nlt(x,t)<[>ly(x,oft)   -  -~[^lx(xfoft)2  +  <$>ly(x,o,t)2)  +  k2 


(26) 
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There  is  a  striking  similarity  between  the  first-order 
and  second-order  problems;  only  the  right  hand  side  of  the 
free  surface  boundary  conditions  differ.   In  the  first-order 
problem  the  free  surface  boundary  conditions  are  homogeneous 
whereas  the  second-order  free  surface  boundary  conditions 
are  non-homogeneous  functions  of  the  first-order  velocity 
potential,  first-order  free  surface  elevation,  and  their 
derivatives . 

The  first-order  potential  function  may  be  represented 
by  a  function  which  is  periodic  and,  therefore,  the  complex 
potential,  u, (x,y) ,  is  defined  as: 


♦x(x,yft)  =  abRe[iu1(x,y)e"lt]  (27) 


where  R  denotes  the  real  part,  a  =  2iTa/L,  L  denotes  the 
wave  length,  and  b  is  an  unknown  real  constant. 


C.   THE  FIRST-ORDER  BOUNDARY  VALUE  PROBLEM 

Since  the  first-order  boundary-value  problem  is  linear, 
the  potential  <j>,  may  be  expressed  as  a  sum: 


♦x  =  <J>*  +  ^  (28) 


I  S 

where  4>,  denotes  the  incident  wave  potential  and  <J>,  denotes 

the  scattering  potential  due  to  the  presence  of  the  cylinder 

In  terms  of  the  time  independent  complex  potentials,  using 

Eqs.  (27-28): 
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ul  =  ul  +  ul  (29) 


Remembering  that  the  incident  wave  potential  represents 
only  the  incoming  wave,  with  all  the  effects  of  the  cylinder 
represented  by  the  scattering  potential,  then  the  incident 
potential  must  satisfy  the  first-order  boundary-value  problem 
when  no  rigid  body  is  present.   Therefore,  <{>..  satisfies 
Eqs.  (17-18)  and  (20-21)  ,  which  represents  simply  the  first- 
order  boundary-value  problem  for  a  train  of  regular  waves. 
The  solution  to  this  well-known  problem  in  terms  of  the 
complex  potential  is: 


I     1  cosh[a(y+h)]   lax  ,^~, 

u,  = ,  )    *  \         e  (30) 

±  a    cosn(an) 


The  relationship  between  v  and  a  is  defined  by 


2- 

v  =  °-JL  =   a  tanh(ah)  (31) 

y 


which,  in  dimensional  form,  using  a  =  ka,  and  k  =  2tt/L  is 


a2  =  gk  tanh (kh)  (32) 


Since  the  solution  for  u..  is  known,  the  boundary-value 
problem  for  u,  may  now  be  established.   Substituting  Eqs. 
(30) ,  (29)  and  (27)  into  the  first-order  problem  given  by 
Eqs.  (17-21),  and  eliminating  n,  between  Eqs.  (20)  and  (21) 
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X 


results  in  a  boundary-value  problem  .for  the  first-order 
scattering  potential  as  follows: 


V2u^(x,y)  =  0  (33) 


u5L(x,-h)  =  0  (34) 

uln(x'Y)  =  cosh  (ah)  [ny  sinhta(y+hH  +  inx  cosh [a (y+h)j  ela 
on  S1(x,y)  =  0  (35) 


u^y(x,0)  -  vu^(x,0)  =  0  (36) 


where  nv  and  n 7   denote  the  x-  and  y-components ,  respectively, 

of  the  unit  normal  vector,  n  =  in   +  in  ,  directed  outward 

'        x    J    y 

into  the  fluid. 


D.   SECOND-ORDER  BOUNDARY-VALUE  PROBLEM 

In  view  of  the  linearity  of  the  second-order  problem, 
the  same  procedure  as  used  in  the  first-order  problem  may 
be  applied  and  leads  to  the  representation  of  the  second-order 
potential  as  the  sum: 


<J>0  =  4>l   +    £  (37) 


where  <f>2  denotes  the  second-order  incident  wave  potential 

5 
and  <f>2  denotes  the  scattering  potential. 
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Considering  the  case  where  there  is  no  body  present, 

S     S 
there  is  no  scattering  wave  and  therefore,  <f>,  =  <f>2  =  0 . 

Additionally,  the  boundary  conditions  on  the  immersed  cylinder, 
Eqs.  (19)  and  (24)  are  not  applicable.   When  the  substitu- 
tions of  (f>.  for  (f>,  and  $~    for  <£~  are  made  into  Eqs.  (22)  , 
(23),  (25)  and  (26),  along  with  Eqs.  (20),  (21),  (27)  and 
(31),  and  upon  eliminating  ru  an<3  rip  between  Eqs.  (25)  and 
(26) ,  the  resulting  boundary-value  problem  for  the  second- 
order  incident  wave  potential  is: 


V2<|)J(x,y,t)  =  0  (38) 


<l>2y(x,-h,t)  =  0  (39) 

(f>2y(x,0,t)  +  v4>2tt(x,0,t)  =  -  |b2(a2-v2)Re[iei2{ax_t)] 

(40) 

The  periodic  solution  to  the  incident  wave  boundary-value 
problem  specified  by  Eqs.  (38-40)  is  known  and  may  be  written 
as: 


4*   =  |ab2vRe[iu2(x,y)e"l2t]         (41) 


where  the  second-order  complex  potential,  u~  ,  is  given  by: 


I,    ,      1  cosh[2a(y+h) ]   i2ax 
u0ix,y;  -  -  ~— -  -, e 


(42) 
sinh  * (ah) 
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Eqs.  (41-42)  are  familiar  in  wave  theory  and  exactly  the 
same  form  as  that  given  by  Ref.  4  for  second-order  Stokes' 
waves. 

Having  developed  a  solution  for  the  second-order  incident 

wave  potential,  <{>2 ,  the  problem  for  second-order  scattering 

IS  IS 

potential  may  be  determined.   <{>,  =  <j),  +  4>,  and  (j>2  =  <t>2  +  $j 

are  substituted  into  Eqs.  (22-26),  and  the  previously  deter- 
mined solutions  for  the  incident  wave  potentials,  Eqs.  (30) 
and  (42) ,  are  utilized.   After  applying  the  relationships  of 
Eqs.  (20),  (21),  (27),  (29)  and  (41)  and  eliminating  n2 
between  Eqs.  (25)  and  (26),  the  resulting  boundary-value 
problem  for  the  second-order  scattering  potential  becomes: 


V2$l(x,ytt)    =   0  (43) 


<J>2    (x,-h,t)    =    0  (44) 

S  3  ab2    Re    r 

<J>,n(xfy,t)    =  -   j  z <nv   cosh[2a(y+h)] 

zn  H    sinh    (ah)     L     x 

-in      sinh[2a(y+h) ])el2(ax~t)J  (45) 

on  S    (x,y)    =   0 

aS     ^s      a2b2  D  r.#l  S   S    .   I  S      c    I   S     fAC. 
*2y  +  *2tt  =  "IT-  Re[l(vulyulyy  +  ululyy  "  6ulyuly    (46) 

1SI  .    I      S  0S2  0S    2v     -i2t, 

+  — u,    u,         -    4u,    u,       -    2u.         -    3u,       )e 
v   ly   lyy  lx   lx  lx  ly 


+  U(x)  on   y  =   0 
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where  U(x)  is  a  non-periodic  function  generated  by  the 
substitution  of  <J>.  into  Eqs.  (25)  and  (26).   For  brevity, 
this  function  is  not  written  out  since  it  is  not  needed  in 

determining  the  second-order  pressures  and  forces. 

S 
The  boundary-value  problem  in  <J>2  given  in  Eqs.  (43-46) 

is  time  dependent  according  to  e     ,  except  for  the  U(x) 

term  in  Eq.  (46).   Therefore,  the  solution  for  <j>2  may  be 

taken  in  the  form: 


<j>2  =  4  ab  vRei[u2(x,y)e_1   +  u2(x,y)]  (47) 


where  the  last  term  denotes  the  time  independent  portion  of 

the  complex  potential.   Separate  boundary-value  problems  for 

S      ~S 

u2  and  u2  arise  from  the  substitution  of  Eq.  (47)  into  Eqs. 

(4  3-4  6) .   However,  as  will  be  demonstrated  in  the  pressure 

S 
and  force  development,  only  the  <j>2   term  is  required.   There- 

fore,  the  time  independent  part  of  <j>2  will  be  zero,  and  as  such 

will  not  be  further  developed;  concentration  will  be  directed 

S 
towards  the  solution  for  u? . 

When  Eq.  (47)  is  substituted  into  Eqs.  (43-46),  the 

resulting  boundary-value  problem  for  the  second-order 

scattering  potential  is: 


V2u2(x,y)  =  0  (48) 


u2y(x,-h)  =  0  (49) 
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s  i         r 

u9    (x,y)    =  2 hv    sinh[2a(y+h)  ] 

zn  sinh4  (ah)     L  y 


+    inx   cosh[2a(y+h) ]]    el2ax      on   S1(x/y)=0        (50) 


U2y(x,0)    -    4vu2(x/0)    =    f*(x)  (51) 


where 


f*(x)    =  §^ulyulyy  +  ululyy  "    6ulyuly  +  ^uIyuIyy  (52) 

2  2 

-    4u,    u,       -    2u,         -    3u,  ,    ]       n 
lx   lx  lx  ly     y=0 


Again,  there  is  similarity  in  form  between  +-he  first-order 
and  second-order  problems;  the  only  difference  in  form  being 
that  the  second-order  free  surface  boundary  condition, 
Eq.  (51),  is  non-homogeneous. 

Since  Eq.  (51)  is  non-homogeneous,  further  use  of  the 
linear  superposition  theory  is  made  by  defining  u~  as  the 
sum: 

u2  =  u2   +  u2  ^53^ 

S°  S* 

where  u~   denotes  the  homogeneous  solution  and  u?   the 

particular  solution  to  the  boundary-value  problem  as  stated 

S*      S° 
in  Eqs.  (48-51).   More  precisely,  u~   and  u2   are  defined  as 

the  solutions  to  the  boundary-value  problems  obtained  from 

the  substitution  of  Eq.  (53)  into  Eqs.  (48-51)  as  follows: 
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s* 
Particular  Solution    (u^    ) 


V2u^*(x,y)    =    0  (54) 


U2*    (x,-h)    =   0  (55) 


U2*(xf0)    -    4vu2*(x,0)    =    f*(x)  (56) 


Homogeneous   Solution 


r,2   S°.         ,         n  (57) 

V    u2    (x,y)    =0 

S° 
u^    (x,-h)    =    0  (58) 


S°  S° 

u*    (x,0)    -    4vu2    (x,0)    =    0  (59) 


,0 


(60) 


vu(x,y)    =  i \n      sinh[2a(y+h)] 

z  sinh    (ah)       L  y 

+    inx  cosh[2a(y+h)]lel2ax   -   u^^y) 

on  S, (x,y)    =   0 


By  the  judicious  division  into  homogeneous  and  particular 

S* 
solutions,  the  non-homogeneous  problem  for  u~   contains  no 

boundary  condition  on  the  cylinder  surface.   This  boundary- 

S* 
value  problem  for  u~   as  stated  in  Eqs.  (54-56)  is  identical 

to  that  associated  with  the  linear  problem  for  the  potential 
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resulting  from  a  harmonic  pressure  variation  of  amplitude 

distribution  f *  (x)  on  the  free  surface  in  water  of  depth  h. 

S° 

The  homogeneous  boundary-value  problem  for  u~  ,  defined  by 

Eqs.  (57-60) ,  is  similar  in  form  to  the  first-order  problem 
given  by  Eqs.  (33-36),  the  significant  difference  being  the 

term  4v  in  place  of  v  which  occurs  in  the  free  surface  boun- 

S° 

dary  condition.   Thus,  the  method  of  solution  for  u?   will 

be  similar  to  that  of  the  first-order  scattering  potential. 


E.   DEFINITION  OF  THE  INCIDENT  WAVE 

Having  developed  the  first-order  and  second-order  boun- 
dary-value problems,  it  is  now  appropriate  to  completely 
define  the  incident  wave  height  and  in  the  process  determine 
expressions  for  the  unknown  constants  b  and  K~ .   Solving 
Eqs.  (21)  and  (26)  for  n,  and  ru '  respectively,  and  then 
substituting  the  results  into  the  expression  for  the  free 
surface  elevation  as  given  by  Eq.  (14)  yields: 


n(x,t)  =  -  E*lt  +  £2[-*2t  +  *lt*lyt  +  *itt*ly       (61> 


l(*lx+  *ly>  +  K2>  +  0(e3» 


where  the  potentials  are  evaluated  at  y  =  0 .   Eq.  (61)  is  an 
expression  of  the  free  surface  elevation  in  terms  of  the 
total  potential,  and  therefore,  includes  the  effects  of  both 
the  incident  and  the  scattering  waves.   Hence,  as  it  is  of 
interest  to  obtain  an  expression  for  the  free  surface  elevation 
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of  the  incident  wave,  the  scattering  potentials  are  set  to 

S     S 
zero,  i.e.  <J>.  =  <J>_  =  0,   Evaluating  Eq.  (61)  using  the  known 

solutions  for  the  incident  wave  potentials,  Eqs .  (30)  and 

(42),  along  with  Eqs.  (27)  and  (41),  then  yields: 

2,2      2 


Jt..  ^    _   ^r,    r„i(ax-t),        „2rV(v~-a*) 


H    (x,t)    =    ebRe[e"vcl"    w]    +    £    [  4v  +   K2  (62) 


+    ab2   cosh  (ah)    [2   +   cosh  (2ah)  ]      R    ,    i2(ax-t),"| 
4  cinh3  r^v^  e  J 


+    0(e3) 


where  r\    (x,t)  denotes  the  dimensionless  surface  elevation 
for  the  incident  wave  with  no  cylinder  present.   If  the  x-axis 
is  placed  at  the  mean  water  line  then  the  constant  term  must 
vanish  and,  therefore, 


k2/  2     2\ 
K2 4^ (6j) 


The  remaining  terms  in  Eq.  (62)  are  time  dependent  periodic 
functions,  the  second-order  term  at  twice  the  frequency  of 
the  first-order. 

Defining  the  dimensionless  wave  height  given  in  Eq.  (7) 
as  the  difference  in  elevation  between  the  crest  and  trough 
of  the  incident  wave,  then  eb  must  represent  the  wave  ampli- 
tude.  The  second-order  term  in  Eq.  (62)  at  twice  the  funda- 
mental frequency  makes  equal  contributions  to  surface  elevation 
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at  both  the  crest  and  trough,  so  that  it  contributes  nothing 
to  the  wave  height.   Thus,  in  terms  of  dimensionless  wave 
height,  b  is  given  by 

H  =  eb  =  H/2a  (6  4) 

where  H  denotes  the  elevation  difference  between  the  wave 
crest  and  trough. 

Expressing  the  incident  wave  profile  in  terms  of  the 
dimensionless  wave  height  yields: 


riI(x,t)  =  H  cos  (ax-t)  (65) 


2 

,  H  a  cosh(ah)  ro  ,     ,  /.-  ^  ,      r~  ,    ,. , 
+  __ [2  +  cosh  (ah)  ]  cos  [2  (ax-t)] 

4  sinh  (ah) 


It  may  be  noted  that  Eq.  (65)  agrees  with  the  expression 
for  the  second-order  Stokes'  wave  given  by,  for  example, 
Ref.  4. 

F.   PRESSURES  AND  FORCES 

The  pressure  may  be  formulated  by  use  of  Bernoulli's 
equation,  arranged  as  follows: 


2      2 

P(x,y,t)  =  -p$7r  -  ~p[$-  +  (J>-  ]  -  pgy  +  pgK  (66) 

^-    ^   x     y 


where  P(x,y,t)  denotes  the  dimensional  pressure  and  p  denotes 
the  fluid  density.   By  use  of  Eqs.  (7),  (13-16),  (27),  (41), 
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(47) ,  and  (64) ,  Eq.  (66)  may  be  reduced  to  dimensionless  form, 
and  carried  to  the  second-order  in  wave  height  as  follows: 


P(x,y,t)  =  -y  -  HaRe[u1e"lt]  (67) 


2     2  r-  2 

i    r6v 

^eL    a  u2         "lx        "ly' 


H    a    L     r6v"  2  2.     -i2t, 
4v~L  -  —r~U°    ~    Ul"   ~   Ul"} 


2    '       2 


lulxl   +  lulyl   +  77  "  !] 

J      a     -■ 


In  Eq.  (67)  p(x,y,t)  denotes  the  dimensionless  pressure 
coefficient  and  is  defined  as: 


p(x,y,t)  =  P(X/^t}  (68) 

pga 


The  first  term  in  Eq.  (67)  represents  the  hydrostatic 
pressure  as  y  is  the  dimensionless  depth  beneath  the  mean 
free  surface  (y  =  0) .   The  second  and  third  terms  are  harmonic, 
the  third  term  having  twice  the  frequency  of  the  second  and 
representing  the  harmonic  second-order  contribution.   The 
remaining  terms  in  Eq.  (67)  are  independent  of  time  and 
provide  the  time-average  or  steady-state  force  components. 

Expressing  the  components  of  the  wave  force  vectors  in 
terms  of  integrals  of  the  pressure  over  the  cylinder  surface 
area,  the  dimensionless  force  coefficients  may  be  written 
as: 


Ci(t)  =  -  yp(x/y/t)ni  dC1    i  =  1,2         (69) 
Cl 


where  the  dimensionless  force  coefficients  in  the  x  and  y 
directions  are  defined,  respectively,  as: 


Pga' 


Cx(t)  =   X_3  (70) 


F  (t) 
C2(t)  =  -£~-  (71) 

pga 


with  F  (t)  and  F  (t)  denoting  the  force  components.   Addi- 
tionally, dC,  denotes  a  dimensionless  differential  arc 
length  along  the  circumference  of  the  cylinder. 
Applying  Eq.  (67)  to  Eq.  (69)  yields: 


C,  (t)  =  -  y    /n.dC.  -  Ha   /R  [u,e  Xt]n.  dCn 
i  r  J    i   1       J    e      1  J  i    1 


(72) 


Cf  Cl 


„2ra2    fl     r6v2        2      2.  -i2t 
-  H  [4^ryRe[-a-u2  ~  ulx  -  uly)e 
ul 


2         2     2 

+  |ulx'   +  |uly!   +  ^   "  1]nidci  +  0(h3) 

a 


As  the  first  term  in  Eq.  (72)  results  from  the  hydrostatic 
pressure  on  the  cylinder,  it  represents  simply  the  buoyancy 
force.   Omitting  the  hydrostatic  force,  the  hydrodynamic 
force  may  be  written  as: 
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C^t)  =   H  F1JL  cos(61±-t)    +  H2  (73) 


2  ^        ? 

+  H  [F_.  cos(<52i-2t)  +  F^.]  +  0 (HJ) 
2i  2i 


where  the  first-order,  second-order,  and  steady-state  force 
coefficients  and  phase  shift  angles  are  defined  by  comparison 
of  Eqs.  (72)  and  (73)  as  follows: 


•   Fiiel  u  ■  %Anidci  (74) 

ci 

T-.      2i    a     //6v        2      2N   ,_  ,~cv 

F2ie     =  4^7  r-/(— u2  "  ulx  "  uly)nidCl  (75) 

F2i  "  £  c  /<  iulx!2  +  lulyi2  +  4  "  «»i*i  "6> 


The  dimensionless  force  coefficients,  F. .  and  F„.  are  real. 

The  first  term  in  Eq.  (73)  represents  the  first-order 
forces  of  the  fundamental  frequency,  a.   The  periodic  portion 
of  the  second  term  in  Eq.  (73)  represents  the  second-order 
contribution  to  the  force  at  twice  the  fundamental  frequency. 
The  last  second-order  term  represents  the  steady-state  or 
time-independent  contribution,  sometimes  called  drift  force. 

Evaluation  of  the  force  coefficients  defined  by  Eqs. 
(74-76)  is  the  main  objective  of  this  thesis.   Therefore, 
it  is  necessary  to  evaluate  the  first-order  and  second-order 
complex  potentials,  u,  and  u~  ,  as  well  as  derivatives  of  u. 
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on  the  surface  of  the  cylinder.   Additionally,  evaluation 
of  u.  and  its  derivatives  on  the  mean  free  surface  (y  =  0) 
is  required. 
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III.   METHOD  OF  SOLUTION  AND  NUMERICAL  PROCEDURES 

A.   SOLUTION  OF  THE  FIRST-ORDER  PROBLEM 

The  use  of  a  Green's  function  to  express  the  solution 
to  the  first-order  boundary-value  problem  was  first  formu- 
lated by  John  [Ref .  5] ,  and  applied  to  submerged  ellipsoids 
by  Garrison  and  Rao  [Ref.  3],   This  method  is  considered  to 
be  applicable,  in  principle,  to  arbitrary  shapes  and  is 
mathematically  the  most  straightforward.   The  Green's  func- 
tion, G,  of  unit  strength  which  satisfies  Eq.  (33)  as  well 
as  the  boundary  conditions,  Eqs.  (34)  and  (36)  is  given  by: 

00 

c(Mf€ff|»v)  =  m  r  -  m  r'  +  2PVo/  [^^n^Sg-!!i^4sr 

■  e^^^(^^>]co.l^ldll  (77) 

"  i  2ah+sinh(2ah)  COsh  [al (y+h)  ]  COsh  [al  (r'+h)  J  COS  [al  1  x"^  3 
where: 


r  =  [(x-o2  +  (y-n)2]2  (78) 


r'  =  [(x-^)2  +  (y+n)2]2  (79) 


The  symbol,  a  ,  is  defined  in  terms  of  h  and  v  as  the 
solution  to  the  equation: 
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F(airh,v)  =  0  (80) 


where  F  is  defined  by: 


F(a  ,h,v)  =  a.  tanh(a  h)  -  v  (81) 


Comparison  of  Eqs .  (31),  (80)  and  (81)  demonstrates  that  a, 
is  clearly  the  equivalent  of  a.   However,  this  will  not  be 
the  case  for  the  second-order  problem,  requiring  the  use  of 
separate  notation.   In  Eq .  (77),  PV  denotes  the  principal 
value  of  the  integral. 

This  Green's  function  was  also  given  in  series  form  by 
John  [Pef.  5]  as: 

2   2 
»  (yk+v  )cos[yk(y+h) ]  -yk|x-£ 

G(x,y;  £,n;v)  =  2 it  2^ 2 — 2 cos[y  (n+h)]e 

k=l  v^^y^l^+v  ) 

ia.  |  x-£  | 
4-rrcosh[a    (y+h)  ]  cosh  [a,  (n+h)  ]e 

-  i  ± i (82) 

2a  h  +  sinh  (2a-,h) 


where  a^  is  as  defined  in  Eq.  (80)  and  the  quantities  y, 
are  defined  as  the  real  positive  roots  of  the  equation: 


K(yk,h,v)  =  0  (83) 


where  the  function  K  is  further  defined  as 
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K(p,h,v)  =  y  tan(yh)  +  v  (84) 


Following  the  Green's  function  method  of  solution,  u, 
is  written  as  the  integral  over  the  cylinder  arc  length, 
C1#  as: 


h  Jh^. 


ui(x,y)  =  2T  yf1(C#n)G(x,y;Cfn;v)  dcx      (85) 
ci 

where  (£,ti)  denotes  points  on  the  immersed  surface,  f. (£,n) 
denotes  the  unknown  source  strength,  and  dC,  =  dC,/a  denotes 
the  differential  arc  length  on  the  cylinder  surface,  made 
dimensionless  with  the  characteristic  length,  a. 

Although  the  solution  to  the  first-order  boundary-value 
problem  as  stated  in  Eqs .  (33-36)  is  given  by  Eq.  (85) ,  the 
source  strength  function,  f.(£,n),  must  be  determined  in 

order  to  evaluate  the  potential.   From  potential  theory, 

g 

the  normal  derivative  of  the  potential,  u, (x,y) ,  on  the 

surface  of  the  cylinder  is: 


uln(x'Y)  =  I   fl(x'Y)  +  2„  r 

Ll 


^   /f1(^,ri)Gn(xfy;?fn;v)dC1  (86) 


where  G  ,  the  normal  derivative  of  G,  may  be  determined  by 
differentiation  of  either  Eq.  (77)  or  (82)  in  a  straight- 
forward manner. 

Applying  the  final  boundary  condition,  Eq.  (35),  leads 
to  an  integral  equation  which  may  be  solved  for  f ..  , 
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I  f^x'rY)    +  -^   /f1(5,Ti)Gn(x,y;C,n,v)  dC1 


(87) 


;osh(ah)  ["ysi"h[a(y+h)1  +  ii^cosh  [a  (y+h)  ]] 


lax 
e 


Th 


e  solution  for  f,  from  Eq .  (87)  may  then  be  used  in  Eq. 

S 
(85)  to  determine  the  potential,  u, . 

B.   SOLUTION  OF  THE  SECOND-ORDER  PROBLEM 

The  similarity  in  form  of  the  boundary-value  problem 
for  the  homogeneous  part  of  the  second-order  potential, 
Eqs.  (57-60),  to  the  first-order  potential  is  now  utilized. 

Since  the  only  significant  difference  occurs  in  Eq.  (59) , 

S° 

where  v  is  replaced  by  4v,  we  may  represent  u„   in  a  form 

similar  to  Eq.  (85)  as: 

S°  1     f 

U2  (X/Y)  =  J^        /f2U,n)G(x,y;£,v;4v)  dC]L  (88) 

Cl 

where  G  (x,y ;  E,  ,n ;  4v)  is  defined  by  replacing  v  by  4v  in  Eqs. 
(77)  and  (82).   Therefore,  a2  is  defined  by: 


F(a2,h,4v)  =  0  (89) 


and  a,  is  replaced  by  a2  in  Eqs.  (77)  and  (82)  also.   In  the 
case  of  Eq.  (82) ,  y,  is  defined  as  the  real  positive  roots  of 
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K(yk,h,4v)  =  0  (90) 


The  source  strength  function,  f2(£,r|),  appearing  in 
Eq .  (88)  is  again  determined  by  applying  the  kinematic 
boundary  condition  on  the  cylinder  surface  as  given  by  Eq. 
(59) .   The  integral  equation  for  f2  may  then  be  given  by: 

|  f2(x,y)  +  JL  yf2(C,Ti)Gn(x/y;^n;4v)  dc1         (91) 


Cl 


h   sinh[2a(y+h) ]    +   in   cosh  [2a (y+h) ] 


i2ax    S*     , 
e     -  u2n(x/Y) 


4 
sinh  (ah) 


on  S.  (x,y)  =  0 


S* 
However,  to  solve  Eq.  (91)  for  f  2  ,  first  a  solution  for  u? 

S  * 

must  be  obtained  and  u~   evaluated  on  the  cylinder  surface. 

2n 

S* 
Thus,  at  this  point  the  solution  to  u2   is  sought. 

S* 
The  complex  potential  u~   is  defined  as  the  solution  to 

the  boundary-value  problem,  Eqs.  (54-56).   The  form  of  this 

problem  is  recognized  as  being  equivalent  to  that  of  fluid 

motion  produced  by  a  periodic  pressure  distribution  (as 

specified  by  f*(E,))    with  frequency  2a   on  the  free  surface.   No 

cylinder  is  considered  to  exist  in  the  fluid  region. 

The  solution  to  this  problem  may  be  formulated  as  an 

integral  over  the  free  surface  extending  from  negative 

infinity  to  positive  infinity.   Using  the  present 
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s* 

dimension less  representation,  u2   becomes: 

U2*(X'Y)  ==  ~h        /f*(OG*(x,y;£,0;4v)  d£  (92) 

C2 

The  Green's  function,  G* ,  for  the  problem  is  given  in  Ref.  8 
Upon  transforming  the  solution  to  the  non-dimensional  form: 

G*(x,yK,0;4v)  =  -2PV   f   cosh  [y  (y+h)  ]  cos  [y  (x-Q  ]dy 

-oo     4vcosh(yh)  -  ysinh(yh) 


4iTcosh  [a9h]  cosh  [a9  (y+h)  ]  cos  [a9  (x-£)  ] 
-i  1 1 1 (93) 

2a2h  +  sinh(2a2h) 


where  a9  is  defined  by: 

F(a2,h,4v)  =  0  (94) 

S* 

The  normal  derivative  of  u9   is  required  in  order  to 

completely  define  the  kinematic  boundary  condition  on  the 
cylinder  as  given  in  Eq.  (60) .   From  Eq .  (92),  this  becomes: 

u2*(x,y)  =  i   /f*(S)Gn(x,y;£,0;4v)  d£      (95) 
C2 

The  derivative  of  G*  in  the  direction  normal  to  the  cylinder 
surface  may  be  obtained  by  differentiation  of  Eq.  (93)  in  a 
straightforward  manner.   Thus  using  the  values  for  f* (£)  , 
as  defined  in  terms  of  the  first-order  solution  only,  as 
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s*     s* 

given  in  Eqs.  (52)  and  (56),  u~   and  u?   may  be  determined. 

S*     S° 

Using  u~   ,  u~   may  then  be  determined  to  complete  the 

second-order  boundary-value  problem  solution. 


C.   NUMERICAL  METHODS 

Determination  of  the  forces  on  the  cylinder  surface 
requires  solving  for  the  potentials  u, ,  u2  and  the  derivatives 
of  u,  at  points  on  the  surface,  as  shown  in  Eq.  (72).   Thus, 
the  problem  is  now  one  of  solving  for  both  the  first-order 
and  the  second-order  scattering  potentials,  which  requires 
solving  for  the  source  strength  functions  f,  and  f~,  as  well 
as  the  free  surface  pressure  distribution  source  strength 
function,  f*.   In  view  of  the  complexity  of  the  equations 

it  is  natural  to  attempt  a  numerical  solution. 

S 

The  first-order  scattering  potential,  u, ,  as  given  by 

Eq.  (85)  is  dependent  upon  the  first-order  source  strength 
function,  f , .   Thus,  it  is  necessary  to  solve  Eq.  (87)  for 
f-,  to  determine  u.  .   A  numerical  solution  may  be  developed 
by  dividing  the  cylinder  surface  into  elements  of  length 
AG  =  2iT/m,  with  the  center  of  each  element  assigned  an  index. 
Since  f, (£,n)  is  recognized  as  a  well-behaved  function  for  a 
smooth  surface  cylinder,  it  is  appropriate  to  define  the 
following: 

a-ii(v)  =  \  /Gn(x.  ,y.  ;C,n;v)  dC,  (96) 

ID       it  Ar,  J    n   i   l  1 


AC, 

id 


i,j  =  1,2,3,  .  . . , m 
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i        cosh 


EW     [VXi'yi'    sinMa(yi+h)]  (97) 


+    in    (x. ,y.)    cosh[a(y.+h) ] 1 

X  1  J-  -L 


lax. 
e        x 


fli   "   VV^i'  (98) 


and  thus  Eq.  (87)  may  be  approximated  by  the  complex  matrix 
equation: 


fli  +  aij(v)flj  =  2hlj      i/j  =  i'2"-"™  <99) 


Upon  evaluation  of  a. . (v)  using  Eq .  (96),  the  inversion  of 

Eq.  (99)  may  be  carried  out  on  the  digital  computer  to 

determine  f. .  at  each  nodal  point  on  the  surface  of  the 
ID 

cylinder. 

For  the  purpose  of  evaluating  a  and  3/  either  Eq .  (77) 

or  Eq.  (82)  may  be  used.   For  evaluation  of  a  and  3  when  i  =  j 

Eq.  (77)  must  be  used  to  allow  separate  treatment  of  the 

logarithmic  singularity  as  r  approaches  0.   The  difficulty 

encountered  with  the  logarithmic  singularity  may  be  overcome 

by  carrying  out  its  integration  analytically.   Garrison 

[Ref.  2]  showed  for  the  calculation  of  a. .  that  the  contribu- 

iD 

tion  of  the  normal  derivative  of  the  In  (r)  in  Eq.  (77) 
integrated  over  the  singular  element  of  A8  is  A6/2,  not  only 
for  the  nodal  point  (x.,y-)  =  (E,,r\)  ,    but  for  all  nodal  points 
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(x.,y.)  on*  the  cylinder  surface.   Additionally,  to  calculate 
3. .  Garrison  developed  a  numerical  approximation  for  the 
integral  of  the  In  (r)  over  the  singular  element  of  AG. 

A  second  singularity  arises  in  the  evaluation  of  the 
infinite  integral  occurring  in  Eq.  (77) .   The  first  term  in 
the  integrand  is  singular  at  p  =  a.   Since  the  numerator 
approaches  zero  like  (y-a)  near  a,  Garrison  [Ref.  2]  sub- 
tracted l/(y-a)  from  the  singular  term  in  the  range  0  <_  \i    <_   2a 
and  added  the  contribution  of  the  integral  of  l/(y-a)  which 
in  this  case  was  zero.   This  converted  the  singular  integrand 
to  a  regular  function.   Applying  this  technique  to  numeri- 
cally evaluate  the  resulting  integral  between  0  and  2a, 
Simpson's  three-eights  rule  is  used  with  an  odd  number  of 
subdivisions,  thus  ensuring  that  the  point  p  =  a  is  not- 
encountered. 

With  f..  determined  from  the  inversion  of  Eq.  (99),  the 
first-order  potential  and  its  derivatives  may  be  evaluated 
on  the  cylinder  surface.   Replacing  the  surface  integrals 
with  summations,  Eq.  (85)  and  its  derivatives  may  be  written 
as: 


Uli  =  Bij(v)flj      i/j  =  1'2'3"-"m  (10°.) 


u?  .  =  6  •  •  (v)f.  .  (101) 

lxi    xi;j     Id 


ulyi  "   Byij(v,flj  (102) 
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s    s 

where  u, • ,  u,  .  ,  etc.  denote  functions  evaluated  at  the 
li'   lyi' 

i   nodal  point  on  the  cylinder  surface.   The  complex 
matrices  in  Eqs.  (100-102)  are  defined  as  follows: 

3ij(v)  =  "2T    jG(xi>Yi''Z,r\;v)aci        i,j  =  1,2, ...,m  (103) 
AClj 


3xij(v)  =  TF  hr    fizl*±'Y±' **'*)&! 


(104) 


AC1D, 


3  .  .  (v)  =  7T-     /G  (x.  ,y .  ;C/n;v)dCn 
yij      2tt  ._  J    y   i  Ji         1 


(105) 


AClj 


The  first-order   incident  potential  and  its  derivatives  may 
be  evaluated  directly  at  each  nodal  point  after  differentia- 
tion of  Eq.  (30)  as  needed,  and  combined  with  the  scattering 
potential  and  its  respective  derivatives. 

To  solve  the  second-order  problem,  the  function  f*(£) 
required  in  Eq.  (95)  and  as  defined  by  Eq .  (52)  must  be 
evaluated  numerically.   Dividing  the  mean  free  surface 
(y  =  0)  in  the  vicinity  of  the  cylinder  into  n  equal  incre- 
ments, f*  is  evaluated  at  each  nodal  point.   The  numerical 
solution  uses  again  Eqs.  (100-105)  with  y.  =  0  for  the 
purpose  of  evaluating  u, .  and  its  derivatives  at  the  mean 
free  surface  nodal  points.   The  incident  first-order  potential 
and  its  derivatives  are  evaluated  directly  at  each  nodal 
point,  and  combined  with  the  scattering  potential  and  its 
derivatives  to  solve  Eq.  (52)  for  f*. 
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Having  determined  f*  at  all  nodal  points  on  the  free 

S* 
surface,  the  boundary-value  problem  for  u?   given  by  Eqs. 

S*       S* 
(54-56)  may  be  solved,  i.e.  u~   and  u^   may  be  evaluated 

on  the  cylinder  by  solution  of  Eqs.  (92)  and  (95)  respec- 
tively.  Expressing  the  integrals  as  numerical  summations  in 
complex  matrix  form  yields: 


u^*,  =  f*a*.(4v)         i  =  1,2, ...,m  (106) 


2ni     j  ij 


j  =  1,2, ...  ,n 


u^*  =  f*e*. (4v)  (107) 


where 


atn(4v)  =  k  P(x.,y.;r/n;4v)  d£ 


(108) 


AC2j 


Hj^    =h   An     /GJ(xi^i'^'4v>  d^ 


(109) 


AC2j 


s*     s* 

Thus,  using  f * ,  both  u~   and  u_   may  be  directly  evaluated 
at  each  cylinder  surface  nodal  point. 

To  solve  the  boundary-value  problem  for  the  second  part 

S° 

of  the  second-order  scattering  potential,  u^  ,  as  specified 


by  Eqs.  (57-60),  Eq.  (88)  must  be  evaluated.   However,  to 

S° 

evaluate  Eq .  (8  8)  for  u~  ,  the  second-order  source  strength 

function,  f  ~  ,  must  be  evaluated  by  numerically  solving 
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the  integral  equation,  Eq .  (91).   Replacing  the  integral 
equation  with  a  complex  matrix  summation: 


f2i  +  f2jaij(4v)  =  2ki       i/j  =  1'2"--'m  <110) 


where  a.-(4v)  is  defined  by  Eq.  (96), 


f2j  =  f2ixy^  (111) 

k.  =  1 [n  (x.,y.)sinh[2a(y.+h)]  (112) 

1    sinh4  (ah)  Li11  *■ 

1  12ax.     s* 
+  inx(xi/yi)cosh[2a(yi+h)  ]Je      -  ^(x..^) 


Once  a.  .  (4v)  is  evaluated,  Eq.  (110)  may  be  inverted  to  obtain 
the  second-order  source  strength  function,  f~,  at  each  nodal 
point  on  the  cylinder  surface.   Stating  Eq.  (88)  in  summation 
form: 


S° 
u2i  =  f2j6ij(4v)        i/j  =  lf2,3'"-'m  (113) 


where  £..(4v)  is  defined  by  Eq.  (103).   Using  Eq.  (113), 

s°    1D 

u2  may  be  determined  at  each  nodal  point  on  the  cylinder 

S* 
surface,  then  combined  with  u.   and  the  second-order  inci- 
dent potential,  u„,  to  evaluate  the  total  second-order 
potential,  u~  . 
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Since  the  first-order  potential,  u1 ,  and  its  derivatives, 

u,   and  u,  ,  as  well  as  the  second-order  potential.  u~,  are 
lx      ly  c  '2' 

now  known  at  each  nodal  point  on  the  surface,  the  first- 
order,  .second-order  and  steady-state  force  coefficients  and 
phase  shift  angles  may  be  determined  using  Eqs .  (74-76).   As 
shown  earlier,  the  integrals  may  be  replaced  by  summations, 
using  nodal  point  values  and  the  arc  length  increment,  A6. 
Once  each  integral  is  evaluated,  then  the  force  coefficient 
becomes  the  absolute  value  or  modulus  of  the  complex  integral 
result  and  the  phase  shift  angle  is  the  angle  whose  tangent 
is  the  imaginary  part  over  the  real  part  of  the  complex 
integral  result. 

D.   COMPUTER  SOLUTION 

In  order  to  utilize  the  method  of  solution  described 
in  Sections  B  and  C ,  a  computer  program  was  developed  to 
carry  out  the  indicated  calculations.   Although  accuracy  was 
a  prime  consideration,  an  equally  important  requirement  was 
to  minimize  computer  run  time.   To  achieve  this,  every  attempt 
was  made  to  utilize  symmetry  and  to  generate  certain  constants 
and  matrices  to  eliminate  redundant  computations. 

The  program  consists  of  a  main  program  that  flows  in  the 
order  of  computation  presented  in  Sections  B  and  C,  along 
with  four  subroutines.   Two  subroutines  to  solve  the  first- 
order  Green's  function,  using  both  forms  of  the  function, 
Eqs.  (77)  and  (82),  were  developed  by  Garrison  [Ref.  2]. 
These  subroutines,  GREEN  and  GREENS  respectively,  have  been 
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modified  to  include  calculation  of  the  first-order  scattering 
potential  derivatives,  evaluation  of  the  first-order  scattering 
potential  and  its  derivatives  on  the  mean  free  surface,  y  =  0, 

and  evaluation  of  the  modified  Green's  function,  G* ,  for 

S*      S  * 
the  determination  of  u0   and  u~   at  each  cylinder  nodal 

2        2n 

point. 

For  elements  of  the  a  and  $  matrices  corresponding  to 
small  values  of  the  parameter  (x~£)  ,  GREEN  is  used  while  for 
larger  values  of  (x-£) ,  GREENS  is  used.   With  the  exception 
of  the  diagonal  elements  in  Eqs.  (96)  and  (103-105)  and  a  few 
surface  nodal  points  in  Eqs.  (108-109),  the  majority  of 
elements  of  a   and  3  are  calculated  by  GREENS.   This  is  most 
fortunate  as  the  series  form  converges  rapidly,  requiring  much 
less  computer  time  than  the  integral  form. 

Subroutine  GEODAT  reads  the  input  geometrical  data, 
generates  the  matrices  h.  and  k.  as  defined  by  Eqs.  (97)  and 
(112)  respectively,  and  calculates  certain  geometrical 
parameters  and  matrices  for  repeated  use  in  GREEN  and  GREENS. 
Subroutine  COMAT  inverts  the  complex  matrix  equations  and 
thus  is  used  to  invert  Eqs.  (99)  and  (110)  to  determine 
f,  and  f~  respectively. 

A  cross-reference  between  text  and  computer  program 
nomenclature  is  given  in  Table  I. 
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TABLE  I:   Computer  Program  —  Text  Symbol  Cross-Reference 


Text 

Computer 
Program 

Text 

Computer 
Program 

a 

A 

n 

X 

ANX(I) 

d 

D 

ny 

ANY  (I) 

fl 

F(I,1) 

ui 

U1(I) 

f2 

PI (1,1) 

ulx 

U1X(I) 

f* 

FS(L) 

uiy 

UIY  (I) 

F 
*11 

Cl(l) 

u,  (surface) 

U1IS 

F 
12 

Cl(2) 

u,   (surface) 
lx 

U1ISX 

F 
21 

C2(l) 

u,   (surface) 

U1ISY 

F 
22 

C2  (2) 

ulvv  ^surface- 

U1ISYY 

FSS 
21 

C3(l) 

s 
u,  (surface) 

U1SS 

pss 

r22 

C3(2) 

u,   (surface) 
lx 

U1SSX 

G 

GIJ,GIJEXT 

u,   (surface) 

U1SSY 

G* 

GIJ,GIJEXT 

5 
ulyy  ^surface) 

U1SSYY 

Gn 

GNIJ,GNJI 

u2 

U2(I) 

G 
X 

GXIJ,GXJI 

X 

X(I) 

Gy 

GYIJ,GYJI 

y 

Y(I) 

G 

yy 

GYY 

a 

ALPHA  (I,  J)- 

h 

H 

B 

BETA (I, J) 

2hi 

HH(I,1) 

"»x 

BETAX(I,J) 

ki 

PK(I,1) 

*y 

BETAY(I,J) 

m 

NPTS 

A6  (cylinder) 

DELTHE 

n 

NSPTS 

AE  (surface) 

DELX 
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Text 


Computer 
Program 


Text 


Computer 
Program 


11 

PHASEl(l) 

y 

12 

PHASEK2) 

uk(v) 

21 

PH£SE2(1) 

yk(4v) 

22 

PHASE 2 (2) 

V 

K 

X(J) 

4v 

n 

Y(J) 

IT 

AMU 
AMU(K) 
AMU4 (K) 

ANU 

ANU4 
PI 
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IV.   DISCUSSION  AND  RESULTS 

A.   SELECTION  OF  COMPUTER  PROGRAM  PARAMETERS 

The  computer  program  was  developed  on  the  IBM-360 
computer  using  FORTRAN  H.   All  subroutines  were  compiled  on 
DATA  CELL  to  minimize  compile  time  and  core  size.   In  order 
to  maximize  accuracy  while  minimizing  computational  time, 
two  key  computer  parameters  were  evaluated  to  determine 
optimum  values;  cylinder  nodal  points  and  free  surface  model 
points. 

Using  the  first-order  solution  only,  the  number  of  nodal 
points  on  the  cylinder  surface  was  varied  from  12  to  60  and 
compared  with  the  results  determined  by  Garrison  [Ref.  2]. 
The  minimum  number-  of  cylinder  surface  nodal  points  that 
provided  accuracy  within  one  percent  of  those  obtained  by 
Garrison  using  60  points,  was  24.   Accordingly,  the  good 
correlation  from  60  points  down  to  24  points  led  to  the 
selection  of  24  nodal  points  on  the  cylinder  surface. 

Since  the  free  surface  integral  was  to  be  carried  out 
from  -°°  to  +°°,  the  next  step  was  to  establish  the  finite 
interval  of  convergence  on  the  surface  as  well  as  the  sub- 
division size.   After  the  outer  limits  of  the  free  surface 
integral  were  determined,  the  subdivision  size  was  varied 
from  4  to  12  8  subdivisions  per  second-order  wave  length 
holding  the  total  interval  constant.   Above  64  subdivisions 
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the  results  did  not  vary  more  than  one  percent,  leading  to 
the  selection  of  64  subdivisions  per  second-order  wave  length 
oh  the  free  surface. 

Additionally,  the  total  interval,  2X,    was  varied  from 
one  to  eight  wave  lengths  in  both  the  positive  and  negative  x 
directions.   The  second-order  force  coefficients  and  their 
respective  phase  shift  angles  converged  to  a  small  value 
varying  periodically.   Convergence  to  this  limit  cycle 
occurred  between  the  first  and  second  wave  length  out  from 
the  center  of  the  cylinder.   Therefore,  a  total  interval  of 
from  -2  to  +2  second-order  wave  lengths  was  chosen  as  adequate 
for  generation  of  the  numerical  results. 

Although  the  values  of  the  second-order  force  coefficients 
and  their  respective  phase  shift  angles  tended  to  converge 
with  increasing  the  limits  of  the  free  surface  integration, 
there  remained,  in  general,  the  small  limit  cycle  as  noted 
above.   The  amplitude  of  the  cycle  was  a  function  of  the 
parameter,   a,  and  was  of  significant  magnitude  only  for 
values  of  a  less  than  0.25.   The  effects  of  a  on  the  limit 
cycle  for  one  of  the  second-order  parameters,  horizontal  force 
coefficient,  is  demonstrated  in  Fig.  (2).   This  represents  a 
plot  of  percent  variation  of  the  force  coefficient  from  the 
mean  versus  the  surface  interval  size,  X/(L/2),  corresponding 
to  values  of  a  of  0.25,  0.20  and  0.15.   Figure  (2)  is 
representative  of  the  behavior  of  all  second-order  force 
coefficients  and  phase  shift  angles  and  demonstrates  the 
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increase  in  the  amplitude  of  the  variation  with  a  decrease 
in  a.   The  amplitude  of  the  periodic  variation  ranges  from 
up  to  fifty  percent  at  a  wave  number  of  0.15  and  depth  of 
d  =  1.4,  to  less  than  two  percent  for  all  values  of  wave 
number  above  0.25. 

B.   RANGE  OF  APPLICABILITY 

In  a  complex  computer  program  of  the  type  developed  to 
carry  out  the  present  numerical  scheme,  certain  limits  with 
respect  to  numerical  stability,  etc.  show  up.   These  computing 
limits  arise  for  various  reasons  such  as  overflows  caused  by 
computing  hyperbolic  functions  of  very  large  numbers, 
numerical  convergence  instabilities,  etc.   While  no  attempt 
was  made  to  investigate  each  of  these  numerical  problems ,  a 
list  of  the  limitations  of  the  particular  program  developed 
in  this  thesis  to  carry  out  the  computations  are  as  follows: 

minimum  a  —  dependent  upon  acceptable  error,  but  0.25 
or  less 

maximum  a  —  deep  water  condition  reached 

minimum  h  —  3 

maximum  h  —  deep  water  condition  reached  for  a  >  0.25, 
h  =  20  for  a  <  0.25 

minimum  d  —  from  1.4  at  h  equal  to  or  less  than  5  down 
to  2.5  at  h  =  20 

minimum  SMIN  —  0.12  when  cylinder  depth  is  other  than 

near  the  free  surface 

maximum  SMIN  —  0.30  when  cylinder  depth  is  near  the 

free  surface 
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C.   RESULTS 

A  representative  set  of  results  for  the  first-order  and 
second-order  force  coefficients  have  been  generated  for  a 
water  depth  of  h  =  5.0,  submergence  depths  of  d  =  1.5,  2.0, 
2.5,  and  3.0,  with  a  ranging  from  0.15  to  1.2.   All  of  the 
numerical  results  presented  herein  were  based  on  24  sub- 
divisions on  the  circular  cylinder  and  64  subdivisions  per 
second-order  wave  length  on  the  free  surface  integration. 
The  surface  integral  was  carried  out  from  x  =  -L  to  +L  since 
this  was  found  to  be  adequate  for  convergence.   Since  the 
second-order  wave  length  is  half  the  first-order  wave  length, 
L,  the  total  interval,  is  equal  to  four  second-order  wave 
lengths. 

The  first-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (3)  and  (4),  respectively.   It  may  be 
noted  that,  in  general,  the  forces  decrease  with  depth  of 
submergence  according  to  expectation  since  the  wave  action 
dies  out  with  depth. 

The  second-order  horizontal  and  vertical  force  coefficients 
are  presented  in  Figs.  (5)  and  (6),  respectively.   Generally, 
the  results  show  the  second-order  effect  to  be  relatively  more 
important  as  the  cylinder  approaches  the  free  surface.   It 
might  be  suspected  from  this  that  the  second-order  contribution 
would  be  much  more  significant  in  the  case  of  surface  piercing 
or  floating  bodies. 
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The  horizontal  and  vertical  steady-state  force  coefficients 
are  shown  in  Figs.  (7)  and  (8) ,  respectively.   These  results 
show  that  the  horizontal  steady-state  force  coefficient  is  very 
small  in  general.   The  horizontal  force  can  be  shown  to  be 
proportional  to  the  momentum  flux  of  the  reflected  wave  and 
the  reflected  wave  is  small.   In  fact,  in  the  case  of  infinite 
water  depth,  Dean  [Ref.  1]  showed  that  the  reflection  was 
exactly  zero  and  accordingly  the  steady-state  horizontal  force 
is  zero. 

It  may  be  noted  that  the  steady-state  vertical  force  is 
positive.   This  results  from  the  fact  that  the  velocities  are 
larger  on  the  top  of  the  cylinder  and  accordingly,  the 
pressures  reduced. 

The  phase  shift  angles  of  the  first-order  and  second-order 
forces  are  shown  in  Figs.  (9)  -  (12) . 

To  demonstrate  the  effect  that  the  second-order  terms  have 
on  the  forces  on  the  cylinder  and  their  respective  phase  shift 
angles,  a  comparison  was  made  between  the  first-  and  second- 
order  results.   Figures  (13)  -  (18)  are  plots  of  the  horizontal 
and  vertical  dimensionless  forces  on  the  cylinder  versus  time 
over  one  complete  wave  cycle  for  both  the  first-order  solution 
alone  and  the  combined  first-order  and  second-order  effects. 
Additionally,  a  plot  of  the  incident  wave  is  included  to 
demonstrate  the  phase  shift  involved  (the  amplitude  of  the 
incident  wave  has  no  significance).   A  mean  water  depth,  h,  of 
5.0,  a  cylinder  depth,  d,  of  1.5  and  a  wave  height,  H,  of 
0.5  were  used,  with  the  wave  number,  a,  assigned  three  values, 
0.25,  0.5  and  1.0. 
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V.   CONCLUSIONS 

A  computer  program  has  been  developed  to  carry  out  the 
numerical  solution  of  the  second-order  wave  —  cylinder 
interaction  problem.   The  first-order,  second-order,  and 
steady-state  force  coefficients  were  determined  for  the 
submerged  horizontal  cylinder. 

An  approximate  method  for  dealing  with  the  second-order, 
nonhomogeneous  free  surface  boundary  condition  was  developed 
which  appears  to  converge  except  at  small  values  of  a,  i.e., 
large  wave  lengths. 
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//  EXEC  FORTHCLGfREGI0N=350K 
//F0RT.5YSIN  DD  * 

THIS  PROGRAM  CALCULATES  WAVE  FORCES  AND  THEIR  PHASE 

SHIFT  ANGLES  FOR  WAVE  INTERACTION  WITH  A  SUBMERGED 

HORIZONTAL  CYLINDER 

A  =  2*PI*CYL  IN3ER  RADIUS/WAVE  LENGTH  =  SIGMA**2* 

CYLINDER  RADIUS/ACCELERATION  OF  GRAVITV 

H  =  MEAN  WATER  DEPTH/CYLINDER  RADIUS 

D  =  CYLINDER  DEPTH/CYLINDER  RADIUS 

NPTS  =  NUMBER  OF  CYLINDER  SURFACE  ELEMENTS  AND  NODAL 

POINTS  FOR  NUMERICAL  EVALUATION 

NSPTS  =  NUMBER  OF  FREE  SURFACE  ELEMENTS  AND  NODAL 

POINTS  FOR  NUMERICAL  EVALUATION 

COMPLEX  GIJ,GNI JtGNJI ,GXIJ,GXJ I , GYIJ, GY J  I , G IJ EXT , GYY , 
1DET,SUM1,SUM2,S'JM3,SUM4,U1I  S,U1ISX,U1ISY,U1ISYY, 
2U1SS,U1SSX,U1SSY,U1SSYY 

COMPLEX  ALPHA (24, 24) , BETA( 24,24) ,BETAX( 24,24) , 
1BETAY(24,24) ,HH( 24,1; ,F(24, 1) ,PK124, 1) ,F1(24,1) , 
2FS(50O)  ,U1(24)  ,  U1X(2<0  ,U1Y(24)  ,U2(24)  ,J<:SC1(24)  , 
3U2SCN1 (24) ,U2SCO( 24) 

DIMENSION  X(24) ,Y(24) ,ANX<24) , ANY ( 24 ) , CHY ( 2 5 ) ,CHY2(25) 
1,SHY(25) , SHY2  125) ,COEFG(200) , CGEFG2(200) , 
2COSAMU( 200,2  5) ,C0SAM2( 200,25) , SINAMU(2J0 ,25), 
3SINAM2( 200,25) , AMU (200) ,AMU4(200) ,SH2Y( 24),CH2Y(24J 

DIMENSION  CH2)  tC2(2)vC3(2)  ,PHASE1(2)  ,  PHASE  2(2) 

COMMON/CPX/HHtPK 

C  CMMON/ V AR / X , Y , A  NX , ANY 

C0MMCN/GSHY/ChY,CHY2,SHY, SHY  2, SH2Y ,CH2Y 

CCMM0N/GMU/CCSAMU,C0SAM2,SINAMU, S INAM2 , AMU , AMU4, COEFG, 
1C0EFG2 

C CMMON/ CON  ST /H, D , DELTHE f SMI N, NPTS, NSPTS 
LrmUN/cr  i /  r  i 

EQUIVALENCE  <Hh( 1 ),F( 1) ) 

EQUIVALENCE  ( PK  (  i  )  ,  F  1  (  1  i ) 

PI=3. 141592 

READ    INPUT    DATA    AND    CALCULATE    REQUIRED    REPEATING 
DATA    AND    ARRAYS 

CALL    GEODAT    ( A, A2 , ANU, ANU4, SH2 AH , SH2AH2 , SHAH, SHAH2 , 
1CHAH,CHAH2, AO, A A, 3B , CC , DD , A02 , A A2 , BB2 , CC 2, CD2 ) 

DO    100     1=1, NPTS 

DO    100    J=i,I 

XV=X( I ) 

YV=Yl I) 

XC=X( J) 

YC=Y( J) 

ANXI=ANX(  I  ) 

ANYI  =  ANY( I  ) 

ANXJ=ANX( J) 

ANYJ  =  ANY( J  ) 

IF(ABS(X( I)-X( J) ) .LT.SMIN)     GO    TO    50 

SHYI=SHY( I ) 

CHYI=CHY( I ) 

SHYJ=SHY(J) 

CHYJ=CHY( J) 

CALL    GREENS     ( A, ANU , S H2 AH, SH AH, CHAH,COSAMU, S INAMU , AMU , 
lCOEFGtSHYI  ,CHYI , SHY  J, CHY J, I, J,XV,YV,XC,  Y C t ANX I , AMY  I , 
2ANXJ, ANYJ,GIJ,GNIJ,GNJI,GXI J,GXJ1,GYIJ,GYJI ,GYY) 

GO    TO    90 

50    CONTINUE 

CALL    GREEN    ( A, ANU, SH2 AH, SHAH,CHAH, AO, AA , BB , CC , DD, 1 , J , 
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1XV,YV,XC,YC,  ANXI  ,  ANYI  ,GIJ,GNIJ,GXIJ,GYIJ,GYY) 

IFU.EQ.J)  GNJI=GNIJ 
IF  (I. EG. J)  GXJI=GXIJ 
IF  (I. EG. J)  GYJI^GYIJ 
IF( I.EQ.J)  GO  TO  90 

CALL  GREEN  ( A , ANU, SH2 AH, SHAH, CHAM, AO, AA , BB , CC , DD , J  ,  I T 
1XC,YC,XV,YV,ANXJ,ANYJ,GIJEXT,GNJI , GX J  1 , GY J  I , GYY ) 

90  CONTINUE 

CALCULATE  THE  FIRST-ORDER  ALPHA  AND  BETA  MATRICES 

ALPHA ( I , J)=( l./PI )*GNIJ*DELTHE 

ALPHA( J, U=< l./PI }*GNJI*DELTHE 

BETA(I,J)=(1./(2.*PI) )*GIJ*DELTHE 

BETA( J,  I)=BETA(I  ,  J) 

BET AX (I , J)=( l./(2.*PI) )*GXI J*DELTHE 

BETAX( J ,  I  )=< l./( 2.*PI ) )*SXJI*DELTHE 

BETAY( I , J)=( l./(2.*PI ) )*GYI J*DELTHE 

BETAY(  J, I)=(1./(2.*PI) )*GYJI*DfcLTHE 
100  CONTINUE 

DO  120  I=1,NPTS 

ALPHA ( I  ,  I) =ALPHA{ I , I ) +CMPLX ( 1 . 0 , 0  .0 ) 

BETAXt I , I )=BETAX( I , I) +  ANX( I >*CMPLX (0.5 , 0 .0) 

BETAYt I , I)=6ETAY( I , I )+ANY(I ) *CMPLX tO .5, 0 .01 
120  CONTINUE 

GENERATION  OF  THE  FIRST-ORDER  DISTRIBUTION  FUNCTION, 
FlI.Di  BY  INVERSION  OF  ALPHA(  I  ,  J  )  *F(  I  ,  1  )  =  HH(I,1) 

CALL  COMAT  ( 24, 1 , ALPHA ,HH , DE T, INDICA) 

THE  HH(Iil)  MATRIX  IS  REPLACED  BY  THE  DISTRIBUTION 
FUNCTION, F( I ,1  ) 

COMPUTATION  OF  THE  FIRST-ORDER  SCATTERING  POTENTIAL 
FUNCTION,  U1SC(I),  AND  ITS  X  AND  Y  PARTIAL  DERIVATIVES 
UISCX(I)  AND  U1SCYU),  AND  COMBINATION  WITH  THE 
INCIDENT  POTENTIAL  COUNTERPARTS  TO  COMPUTE  THE  TOTAL 
POTENTIAL  FUNCTION  AND  ITS  X  AND  Y  DERIVATIVES 

DO  155  1=1,NPTS 

SUM1=(0. 0,0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

DO  150  J=1,NPTS 

SUM1=SUM1+F( Jf 1) *BETA< I , J) 

SUM2=SUM2*F( J,1)*BETAX(  I  ,  J) 

SUM3=SUM3+F<  J, 1 ) *BETAY( I , J) 
150    CONTINUE 

UK  I  )=SUM1-(CHV(  I  )/(  A*CHAH)  )  *CE  X  P(  CMPLX  (  0  .  0  ,  A*X  (  I  )  )  ) 

U1X(I)=SUM2-CMPLX( 0 . 0, 1 .0 ) *CHY ( I ) *CEXP(  CMPLX( 0.0, 
1A*X(  I  )  )  J/CHAH 

U1Y( I)  =  SUM3-SHY( I ) #C EX P( CMPLXt 0. 0 , A*X  (  I  )  ) )/CHAH 
155  CONTINUE 

EVALUATION  OF  THE  FPEE  SURFACE  PRESSURE  DISTRIBUTION 
SOURCE  STRENGTH  FUNCTION,  FS(L) 

DELX=0.015625*PI/A 

WRITE  (6,180)  DELX 
180  FORMAT  ( 5X// 5X , • DE LX= • , F 1 0. 5/// ) 
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SP=0.0 
CST=.5*DELTHE/PI 

IC0UNT=1 

DO    250    L=1,NSPTS 

SUM1=(0. 0,0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

SUM<t=(0  .0,0.0) 

00    245     I=1,NPTS 

XV=SP 

YV=0.0 

XC=X( I) 

VC=Y( I ) 

ANXI=0.0 

ANYI=1.0 

ANXJ=ANX(  I  ) 

ANYJ=ANY( I ) 

SHYI=SHY(25) 

CHYI=CHY( 25) 

SHYJ=SHY( I ) 

CHYJ=CHY( I ) 

IF     (ABS(XV-XC) .LE.SMIN)    GO    TO    240 

CALL  GREENS  ( A , ANU , Sn2 AH, SHAH, CHAH,COSAMU, S INAMU , AMU, 
ICOEFGf  S'HYI  ,  CHYI,  SHY  J  ,  CHY  J  ,  2  5  ,  I  ,  XV  ,  YV  ,  XC  ,  YC  ,  AN XI  ,  ANYI  , 
2ANXJ,ANYJ,GIJ,GNIJ,GNJI , GX IJ , OX J  I  ,GYI J , GY J I ,GYY) 

GO    TO    241 

240  CONTINUE 

CALL    GREEN    ( A , ANU , SH2 AH , SHAH,CHAH, AO, AA , BB , CC , DO , 2 5, I , 
1XV,YV,XC,YC,  ANXI  ,  ANY  I  ,G  I  J  ,  GNU  ,  GX  I  J ,  GY  IJ  ,  GY  Y ) 

241  CONTINUE 

SUM1  =  SUM+CST*SIJ*F(  1,1) 

Si  mi  ^ r*ii»to./~r'"r-.c'-*w    t     i«ur  f    t         i    % 
UI'i^-OUl'liToo   |    'r^A   i  J''T   t    1   f   l) 

SUM3=SUM3+CST*GYIJ*F( 1,1) 

SUM4=SUM4+CST*GYY*F<  1,1) 
245    CONTINUE 

U1IS=(-1./A)*CEXP(CMPLXO.O,A*SP)  ) 

U1ISX=CMPLX(0.0,-1.0)*CEXP(CMPLX(G.O,A*SP) ) 

UlISY=-TANH<A*H)*CEXP(CMPLX<O.Of A*SP) ) 

U1ISYY=-A*CEXP(CMPLX(0.0, A*SP) ) 

U1SS=SUM1 

U1SSX=SUM2 

U1SSY=SUM3 

U1SSYY=SUM4 

FS( L)  =  ( 2 . *A/ ( 3 . *  ANU ) ) *( Ul I S*U1 SSYY+U1 1 S YY*U1SS+U1 SS* 
1U1SSYY-6.*U1ISY*U1SSY-3.*U1SSY*U1SSY-4.*UIISX*U1SSX 
2-2.*UiSSX*UlSSX) 

IF    ( ICOUNT.EO.O)     GO    TO    248 

SP=FL0AT(L+l)*DELX/2. 

ICOUNT=0 

GO    TO    250 
248    SP=-SP 

IC0UNT=1 
250    CONTINUE 

CALCULATION    OF    THE    PARTICULAR    SOLUTION    PORTION    OF     THE 
SECOND-ORDER    SCATTERING    POTENTIAL    AND    ITS    NORMAL 
DERIVATIVE,    U2SCKI)     AND    U2SCN1U) 

CST=.5*DELX/PI 

DO    350     1=1 ,NPTS 

SUM1=( 0.0,0.0) 

SUM2=(0. 0,0.0) 

SP=0.0 

IC0UNT=1 

DO    345    L=1,NSPTS 

XV=X(I) 
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YV=Y(I ) 

XS  =  SP 

YS=0.0 

ANXI  =  ANX(  I  ) 

ANYI=ANY(I  ) 

ANXJ=0.0 

ANY J= 1.0 

SHYI=SHY< I ) 

CHYI=CHY( I ) 

SHYJ=SHY<25) 

CHYJ=CHY(25) 

IF  (ABSUV-XS) 


LE.SMIN)  GO  TO  340 


CALL  GREENS  ( A2 , ANU4 , SH2AH2 » SHAH2 ,CHAH2 , C0SAM2, SI NAM2, 
1AMU4,C0EFG2,SHYI , CHY I , SHY  J , CHY J , I , 25 t XV , YV , XS , YS , 
2ANXI vANYIf ANXJ , ANY J , G IJ , GNI J , GN J  I , GX I  J, GX J  I ,GYI J , GYJ  I , 
3GYY) 

GO  TO  341 
340    CONTINUE 

CALL  GREEN  ( A2 , ANU4, SH2AH2 , SHAH2 , CHAH2, A02, AA2 , BB2 , CC2 


l»DD2,I t 
2GYY) 


25,XV,YV,XS,YS,ANXI  »  ANY  I  »  GIJ  ,  GNU  ,  GX  IJ  ,  GY  IJ  , 


341    CONTINUE 

SUM1=SUM1+CST*GI J*FS<  L) 
SUM2=SUM2+CST*GNI J*FS(L) 
IF    (  ICOUNIT.F-Q.O)     GO    TO    290 
SP=FL0AT(L+l)*DELX/2. 
ICOUNT=0 
GO    TO    345 

290    SP=-SP 

IC0UNT=1 

345    CONTINUE 

inrri   /   i   \  -ciimi 
utoi/i  i  l  i  -JUni 

U2SCNK  D-SUM2 
350  CONTINUE 


CALCULATION  OF  THE  HOMOGENEOUS  SOLUTION  PORTION 
OF  THE  SECOND-ORDER  SCATTERING  POTENTIALt  U2SC0( I ) 

DO  500  1=1 ,NPTS 
DO  500  J=l ,1 
XV=X( I) 
YV=Y(I) 
XC=X( J) 
YC=Y( J) 
ANXI=ANXl I ) 
ANYI=ANY( I ) 
ANXJ=ANX( J ) 
ANYJ=ANY( J ) 
IF(ABS(X( I)-X(  J)  ) 
SHYI=SHY( I ) 
CHYI=CHY( I ) 
SHYJ=SHY( J) 
CHYJ=CHY( J ) 

CALL  GREENS  ( A2 , ANU4 , SH2AH2 , SHAH2 , CHAH2 , COS AM2, S I NAM2, 
1AMU4,C0EFG2,SHYI , CHY I , SHY  J , CHY J t I t J » XV t Y V , XC , YC , ANX I , 
2  ANY  I t ANXJ, ANY J i GI J, GNI J,GNJ I , GX  I  J  ,  GXJ I , GY I  J , GY J [ , GYY ) 


LT.SMIN)    GO    TO    450 


GO  TO  490 
450  CONTINUE 

CALL  GREEN  (A2, 
1 ,DD2, It JtXVtYV, 
2GYY) 


ANU4.SH2AH2, SHAH 2 , CHAH2 t A  02 , AA 2 , BB 2 ,CC2 
XC,YC , ANX I ,ANYI , GIJ , GNI J , GX I  J, GYI J , 


IF(I.EQ.J)  GNJI=GNIJ 
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IF( I.EQ.J)  GO  TO  490  . 

CALL  GREEN  ( A2 , ANU4, SH2AH2 , SHAH2 , CHAH2 , A02 , AA 2 , BB 2  ,CC2 
1,DD2, J,  I  ,XC,YC,XV,YV,ANXJ,ANYJ,GIJEXT,GNJI ,GXJI,GYJI , 
2GYY) 

490  CONTINUE 

ALPHA ( I  ,  J)=( l./PI )*GYIJ*DELTHE 

ALPHA (J , I)  =  ( l./P  I )*GYJI*DELTHE 

BETAdt  Jh=(l./(2.*PI)  )*GIJ*DELTHE 

BETA( J, 1)=BETA( 1 , J ) 
500  CONTINUE 

DO    510     l=l,NPTS 

PK(I,1)=2*(PK(I,1)-U2SCN1(I  )) 
510    CONTINUE 

DO    520     1=1 ,NPTS 

ALPHA ( I  ,  I)=ALPHA(  I , I)+CMPLX( 1.0,0.0) 
520  CONTINUE 

CALL  COMAT  ( 24, 1 , ALPHA, PK , DET, I NDICA ) 

DO  540  I=1,NPTS 

SUM=(0. 0,0.0) 

DO  530  J=1,NPTS 

SUM=SUM  +  FIU,1  )*BETA{ I, J) 
530  CONTINUE 

U2SC0( I )=SUM 
540  CONTINUE 

EVALUATION  OF  THE  TOTAL  SECOND-ORDER  POTENTIAL,  U2(I) 

DO  550  1=1 ,NPTS 


4)  )  )* 


UU  3PU   1  =  I  ,  PJ  K  1  O 

U2(I)=U2SC1( I)+U2SC0( I )-(CHY(I ) / ( 2.*A* ( SHAH** 
1CEXP(CMPLX(0.0,2.*A*X( I ) ) ) 
550  CONTINUE 

WRITE  (6,560) 
560  FORMAT  ( 5X//9X , • I •  ,  ] 5X , • Ul ( I )  ' , 24X , ' U1X ( I ) •  ,25X f 
1'UIYU ) ' ,24X, »U2( I ) *//) 

DO  580  I=1,NPTS 

WRITE  (6,570)  I, UK  I  )  ,U1X(  I  )  ,U1Y(I),U2(  I  ) 
570  FORMAT  ( 5X , I  5 , 4( F 1 6 .6 , F 14.6 ) ) 
580  CONTINUE 

EVALUATION  OF  THE  FIRST-ORDER,  SECOND-ORDER  PERIODIC, 
AND  SECOND-ORDER  STEADY  STATE  FORCE  COEFFICIENTS  AND 
THE  PERIODIC  FORCE  PHASE  SHIFT  ANGLES 

SUM1=(0. 0,0.0) 

SUM2=(0. 0,0.0) 

SUM3=(0. 0,0.0) 

SUM4=(0.0,0.0) 

SUM5=0.0 

SUM6=0.0 

DO  720  I=1,NPTS 

SUM1=SUM1+U1 ( I)*ANX( I )*DFLTHE 

SUM2=SUM2+U1 (I)*ANY(I)*DELTHE 

SUM3  =  SUM3+(  (6.*ANU*ANU*U2<  1 )  /A  )  -U1X(  I  )  *U1  X  (  I)-U1Y(  I  J 
1*U1Y( I ) )*ANX( I)*DELTHE 

SUM4=SUM4+( (6.*ANU*ANU*U2( I ) /A ) -U1X( I ) *U IX ( I ) -Ul Y ( I ) 
1*U1Y(  I )  )*ANY( I)*DELTHE 

SUM5=SUM5+(CABS(U1X(  I ) *U1 X ( I ) ) +CABS (U 1Y (  I  )*U1Y(I)  ) 
1+ANU*ANU/A**2-1 . 0)*ANX( I)*UELTHE 

SUM6=SUM6+(CABS(U1X(  I)*U1X(  I  )  )  +CABS(U1Y(  I  )*U1Y(I  )  ) 
1  +  ANU*ANU/A**2-1  . 0 ) *ANY( I ) *D EL THE 
720  CONTINUE 

CI  (1)  =  CABS(SUM1*A) 

C1(2)=CABS(SUM2*A) 

C2(1)=CABS(SUM3*A*A/(4.*ANU) ) 
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C2(2)=CABS(SUM4*A*A/(4.*ANU) ) 

C3(  1)=SUM5*A*A/(4.*ANU) 
C3(2)=SUM6*A*A/(4.*ANU) 

PHASE 1( 1)=ATAN2(AIMAG(SUM1*A),REAL(SUM1*A)) 

PHASE1(2)=ATAN2(AIMAG( SUM2* A ) , RE AL ( SUM2* A ) ) 

PHASE  2 (  1)=ATAN2( A  I  MAG ( SUM3* A*A/ ( 4. *ANU )  )  , RE AL ( SUM3* A*A 
1/(4.* ANU)  )  ) 

PHASE2  1 2)=ATAN2( A  I MAG( SUM4* A*A/ ( 4 .*ANU) ) ,  RE AL( SUM4*A*A 
1/(4. *ANU) ) ) 

WRITE  (6,730)  C 1 (  1  ) , C 1( 2) , C 2 { 1 )  , C2 ( 2) , C 3 ( 1 )  ,C3( 2 ) , 
1  PHASE  1 ( 1 )  , PHASE  1 (2),PHASE2(1),PHASE2(2) 
730  FORMAT  (5X///5X, • CI ( 1 ) =' , F8 .5, 5X , ■ Ci ( 2) = • , F8.5//5X,  ■ 
l'C2(l)=',F8.5,5X,«C2(2)=,,F8.5//5X,,C3(l)=,,F8.5,5X, 
2,C3(2)=,,F3.5//5X,,PHASE1(1)=,,F8.5,5X,,  PHASE  1(2)  =  ' , 
3F8.5//5X, • PHASE  2 ( 1  )=•  ,F8.5,5X,  •PHASE2(2)=»  , F8.5//) 

STOP 

END 

THIS    SUBROUTINE    READS    THE    INPUT    GEOMETRICAL    DATA    AND 
CALCULATES    THE    FIRS1     200    ROOTS     OF    AMU*) AN ( AMU*H)     + 
ANU    =    0    AMD    AMU4*TAN( AMU4*H)     +     ANU4    =0.        IT    ALSO 
GENERATES     ARRAYS    OF    CERTAIN    FUNCTIONS    AND    COEFFICIENTS 
USED    IN    GREEN    AND    GREENS    AS    WELL    AS    THE    HH    AND    PK 
MATRICES 

SUBROUTINE    GEODAT     ( A , A2 T ANU , ANU4, SH2AH , SH2AH2 » SHAH, 
lSHAH2,CHAh,CHAH2, AO, AA,bB,CC,DD,A02,AA2 , Bb2 »CC2f DD2 ) 

COMPLEX  HH(24, 1) ,PK(24,1) 

DIMENSION  X<24) ,Y( 24) , ANX( 24) , ANY< 24)  ,CHY(25) ,CHY2(25) 
If SHY(25) , SHY2(25) , CC EFG ( 200 ) , COEFG2 ( 200 ) , AMU (200) , 
2AMU4(20  0) , COSAMJ(20C,25) ,COSAM2( 200,25) , SINAMU( 200,25) 
3,SINAM2(2  00,25) ,SH2Y( 24) ,CH2Y(24)  ,XX(24) 

CGMMON/CPX/HH,PK 

CCMMON/VAR/X,Yt ANX.ANY 

C  CMMON/ CS  H YA,  HY , C  HY  2 , SHY ,SHY2,SH2Y,CH2Y 

COMM0N/GMU/CDSAMU,CUSAM2f SINAMUf SINAM2f AMU , AMU4 , CGEFG, 
1C0EFG2 

COMMON/CON ST/H,D,DELTHE,SM I N,NPTS,NSPTS 

COMMQN/CPI/PI 

RE AD (5, 5)     A,H,D, SM IN, NPTS , NSPTS 
5    FORMAT     (4F10. 5,214) 

ANU=A*TANH(A*H) 

SH2AH=SINH(2.*A*H) 

CHAH=COSH(A*H) 

SHAH=SINH( A*H) 

AO=TANH(A*H) 

AA=1./CHAH**2 

BB=-H*SKAH/(CHAH**2) 

CC=-H*H*( i.-SHAH**2)/(3.*CHAH**3) 

DD  =  H**3=iSHAH*(  2  .  *CHAH**  2  +  3  .  *  (  1.- SHAH**  2)  )  /  <  CHAH**4*  12  ) 

DELTHE=2.*PI/NPTS 

THETA=DELTHE/2. 

DO    15    1=1, NPTS 

X(I)=COS(THFTA) 

Y( I)=SIN(THETA)-D 

ANX( I)=X(  I  ) 

ANY( I )=Y(  I  )+D 

SHY(  I  )=SINH(  A*(H+Y{  I)  )  ) 

CHY(  I  )  =  COSh( A*(H  +  Y( I )  )  ) 

SH2Y(  I)=SINH(2.*A*(H  +  Y(  I  )  )  ) 

CH2Y( I )=COSH(2.*A*(h+Y( I) ) ) 

HH< I,1)=(2./CHAH)*CNPLX(ANY(  I ) *SHY ( I )  , A  NX (  I )*CHY(  I)  )* 
1CEXP(CMPLX(0.0,     A*X(I))) 

PK(  1,1) =( ANY(  I)*SH2Y( I  )  *CMPLX( 0  .0 , 1. 0) * ANX ( I)*CH2Y( I ) ) 
1*CEXP(CMPLX(0.0,2.0*A*X(1 ) ) )/SHAH**4 

THETA=THETA+DELTHE 
15    CONTINUE 

SHY(25) =SHAH 

CHY(25)=CHAH 
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B=ANU*H 
DO    50    K=i,200 
XX(1)=PI*K 
DO    25    1=1,20 
I  1  =  I 

YY=XX( I ) 

XX(H-1)*XX(1)-ATAN2(3,YY) 

IF(A8S(  (XX( I+lJ-XXl I  J )/XX( I +1) )  .LT..0031 )    GO    TO    30 
25    CONTINUE 

IF    ( II.&E.20)    GO    TO    27 
GO    TO    30 

27  WRITE    (6,28) 

28  FORMAT     (5X,'AMU    ROOT    DOES    NOT    CONVERGE'//) 
30    CONTINUE 

AMU(K) =XX(I I )/H 

C0EFG(K)=2.*PI*(AMU(K)*AMU< K) + ANU*ANU) /  ( ANU*AMU(  K) -H5 
1AMU(K)*< AMU(K)*AMU(KJ+ANU*ANU) ) 

NPTS1=NPTS+1 

DO    40    I=1,NPTS1 

IF    ( I .EC.NPTS1 )     GO    TO    35 

COSAMIM  K,I  )=COS(  AMU(K)*(H«-Y(  I  )  )  ) 

SINAMU(K,  I)=SIM(AMU(;<)*(H+Y(I)  )  ) 

GO    TO    40 
35    CCJSAMUtK,  I)=CGS(AMU(K)*H) 

SINAMU(  K,lJ=SIN(AMU(K)*h) 
40  CONTINUE 
50  CONTINUE 

ANU4=4.*AMU 

B2=ANU4*H 

DO  70  K=l,200 

XX(1 )=PI*K 

DO  55  1=1,20 

I  I  =  I 

YY=XX(  I  ) 

XX ( 1  +  1) =XX(1)-ATAN2(82,YY) 

IF(ABS( <XX< I+1)-XX(I> )/XX{  I+1J ) .LT..0001 J  GO  TO  60 
55  CONTINUE 

IF  (  II.GE.20)  GO  TO  57 

GO  TO  6C 

57  WRITE  (6,58) 

58  FORMAT  (5X,«AMJ4  ROOT  DOES  NOT  CONVERGE1//) 
60  CONTINUE 

AMU4(K)=XX< I  I )/H 

C0EFG2(K)=2.*PI*(  AMU4(iU*AMU4(K)  +ANU4*ANU4)  /(  ANU4* 
1AMU4(K) -H*AMU4( K ) *( AMU4 ( K ) *AMU4 < K) +ANU4*ANU4) ) 

NPTS1=NPTS+1 

DO  65  I=1,NPTS1 

IF  (I.EG.NPTS1)  GO  TO  68 

C0SAM2(K,I )=COS( uMU4 ( K ) - ( H+Y( I ) ) ) 

SINAM2 (K, I)=SIN( AMU4 ( K ) * < H+ Y( I  )  )  ) 

GO  TO  65 
68  C0SAM2( K,I )=COS( AMU4( K)-H) 

SINAM2C K,  I)=SIN( AMU4(K)*H) 
65  CONTINUE 
70  CONTINUE 

XX( 1 )=ANU4 

DO  80  1=1,20 

I  I  =  1 

Y2=XX( I ) 

XX(I+1)=4.*ANU/TANH(Y2*H) 

IF  (A3S( lXX(  I+1)-XX( I ))/XX{  1  +  1) )  .LT. 0.3001)  GO  TO  85 
80  CONTINUE 

WRITE  (6,32) 
82  FORMAT  (5X,'A2  ROOT  DOES  NOT  CONVERGE'//) 
85  CONTINUE 

A2=XX( I  I) 

SH2AH2=SINH(2.*A2*H) 

CHAH2=COSH(A2*H) 

SHAH2=S INH( A2*H) 

A02=TANH(A2*H) 

AA2=1./CHAH2**2 

BB2=-H*SHAH2/(CrtAH2**2) 
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CC2=-H*H*( l.-SHAH2**2)/(3.*CHAH2**3) 
DD2=(H**3)*SHAH2*( 2. *CHAH 2** 2+3 . *( l.-SHAH2**2))/ 
1 ( (CHAH2**4)*12. ) 
DO  90  1=1 ,NPTS 
SHY2( I)=SINH(A2*(H+Y( I))) 
CHY2(I)=COSHU2*(H+Y<  I)  )) 
90  CONTINUE 

WRITE     (6,95)     A,A2 ,D,H  ,ANU,ANU4,NPTS,SMIN,NSPTS 
95    FORMAT     (5X///3X, ' A=« , F10.5,4X, • A2=« ,F1D.5,4X, 'D=« , 

, F10.5,4X,  '  A;vJ  =  «  ,F10.5,4X,«  ANU4=«  ,F10.5/ 
, I3,4X, ' SMIN=«  , F10.5,4X, •NSPTS  =  i ,  13///) 
RETURN 
END 

THIS    SUBROUTINE    CALCULATES    G    AND    ITS    DERIVATIVES, 
GXtGY,GN,     AND    3VY,     3Y    USE    OF    THE     INTEGRAL    FORM    FOR    THE 
CASE    (XII)    -    X(J))     LESS    THAN    SMIN 

SUBROUTINE    GREEiM     (  A  ,  ANJ  ,  SH2  AH,  SH  AH.CHAH  ,  AO,  AA  ,BB  ,  CC  ,DD 
1  ,  I  ,  J,X, Y,XI ,ETA, ANX,ANY,G,GN,GX,GY,GYY) 

COMPLEX    G,GN 

COMPLEX    GX,GY,GYY 

DIMENSION    TEST ( 200) ,TESTT( 100)  , TESTTT( 100)  ,SUMOX( 15)  , 
lSUMOYi 1 5) ,N\N(15 ),TESTYY(200) 

CCMMON/CONST/H,0,DELTHE,SMIN,NPTS,NbPTS 

COMMON/CPI/PI 

PI (Y,ETA,XX,AMU)=COSH(AMU*(Y+H) )*COSH(AMU*( ETA+H) )* 
1C0S(AMU*XX)/(C0SH(  AMU*H)**2) 

P2X(Y,ETA,X,XI,AMJ)=-AMU*C0SH( AMU*(ETA+H) )*SIN( AMU* 
1 (X-XI ) )*COSH( AMU* (Y+H) )/(COSH( AMU*H)**2) 

P2Y(Y,ETA,X,XI , AMU) =AMU*COSH( AMU* (ETA+H)  )*S  INH< AMU* 
1 (Y+H) ) *COS(AMU*(X-XI ) ) / ( C OSH( AMU*H )**2 ) 

1( ETA+H) )*COS(AMU*(X-XI) ) /{ C0SH( AMU*H)**2 ) 

QHYiETA,  XX,  AMJ)=-P1  (        Y  ,  ET  A,  XX  ,  AMU)  *  (  AMU-A)/(AMU* 
ITANH(AMU*H)-ANU) 

Q2X(Y,ETA,X,XI,AMU)=-P2X(Y,ETA,X,XI,AMU)*(AMU-A)/(AMU* 
1TANH(  AMU*H)-A\IJ) 

Q  2Y ( Y , ETA , X , XI , AMU ) =-P  2 Y ( Y , ETA , X , X I , AMU ) * ( A  MU-A) / ( AMU* 
1TANH( AMU*H)-ANU) 

Q2YY( Y, ETA,X,XI, AMU)=-P2YY (Y,ETA,X,X1 , AMU ) * ( AMU- A ) / 
1( AMU*TANH( AMU*H) -ANU) 

Q10(Y ,ETA,XX)=-COSH(A*(Y+H) )*C0SH( A*(ETA+H)  )*C0S( A*XX) 
1*A/(C0SH( A*H)**2*< ( A-A-ANU*ANU)*H+ANU) ) 

Q20X(Y,ETA,X,XI ) = A*C 3SH ( A* ( ET A  +  H ) ) *S I N( A* ( X-X I ) ) *C0SH 
1 ( A*( Y+H ) ) *A/(C0SH( A*H)**2*( ( A*A-ANU*ANJ ) *H+ANU) ) 

Q20Y( Y,ETA,X,XI ) =-A*C0SH( A* ( ETA+H) ) *S I NH ( A* ( Y+H ) ) *C0S 
1 ( A* (X-X J ) )*A/ (C0SH(A*H)**2*( ( A* A- ANU* ANU ) *H+ANU) ) 

Q20YY(Y,ETA,X,XI ) =-A*A*CQSH ( A* ( Y+H) ) *C3S H( A* ( ET A+H) ) * 
1C0S(A*(X-XI))*A/(C0SH(A*H)**2*( ( A* A- ANU* ANU ) *H+ANU ) ) 

Q1S(Y,ETA,XX,A'1J)=-Pi  (  Y  ,  E  T  A  ,  XX  ,  AMU)  /  (  A0  + AMU*H*  (  AA  +  BtJ* 
1 (AMU- A) +CC*( AMU-A)**2+DD*( AMU- A )**3) ) 

Q2XS(Y, ETA,X,XI, AMU ) =-P2X ( Y , ET A , X ,X I , AMU )/( A0+AMU*H* 
KaA+BB*  (AMU- A)  +CC-(  AMU- A ) **2+DD* ( AMU-A ) **3)  ) 

Q2YS( Y, ETA,X,XI, AMU ) =-P 2Y ( Y , ETA , X , X  I yAMU)/( A0  +  AMU*H* 
1 (AA+BB* ( AMU-A)+CC*( AMU-A)**2+DD*( AMU- A)** 3) ) 

Q2YYS( Y,ETA,X,XI  , AMU)=-P2YY ( Y , ET A , X , X  I , AMU ) / ( A0+ AMU*H* 
KAA  +  BB*  (  AMU-A  )  +:  C*(  A  >1U-A  )  **  2  +  DD*  (  AMU- A  )**3)  ) 

FUN1 ( Y , ET  A , X  X , AMU ) =EX  P ( -A  MU  *H ) *  S I NH  ( AMU*  ET A ) *SINH 
1( AMU*Y)*COS( AMU*XX) /( AMU*C0SH( AMU*H ) ) 

FUN2( Y, ETA,XX)=4.*3.14159*C0SH(A*(Y+H))*C0SH(A* 
1( ETA+H) )*C0S(A*XX)/(2.*A*H+SH2AH) 

FUNR(X, Y, XI, ETA, ANX,ANY)=( (X-X I ) *ANX+ ( Y +ET A ) *ANY ) / 
1( (X-XI)**2+(Y+ETA)**2) 

FUN3X( Y,ETA,X,XI , AMU ) =EXP ( - AMU*H ) *S INH( AMU* ET A) *S INH 
1  (  AMU*Y  )  *  S  I  iM(  AMU*  (  X-X  I  )  )  /COSH  (  AMU*H  ) 

F  UN 3 Y  (  Y  ,  E T A  ,  X  ,  X  I  ,  A  M U )  =- :"  X  P  (  - AMU*H  )  * S  I  NH  (  AMU*E TA  )  *C OS H 
1 ( AMU*Y)*COS< AMU* (X-X I ) ) /C OS H( AMU*H) 

FUN3YY ( Y , E  TA , X , X  I , AM J ) =-EX  P {- AMU*H) * S I NH ( AMU*  ETA ) * 
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1SINH<  AMU*Y)*AMU*AMU*COS< AMU*<X-XI  ).)  /COSH  (  AMU*H) 
FUN4(Y,ETA,XX,ANX,ANY,SIGN)=4.*3.  141 59* A*COSH ( A 
1 ( ETA+H) )*( ANY*SINH( A*( Y+H) ) *C3S( A* XX) -ANX*COSH( 
2(Y+H))*SIN(A*XX)*SIGN)/(2.*A*H+SH2AH) 


EVALUATION  OF  THE  FINITE  INTEGRAL  IN  G  TO  DETERMINE 
THE  SIZE  OF  THE  SUBDIVISIONS 

XX=ABS(X-XI ) 

IF(X.LT.XI)  SIGN =-1.0 

IF(X.GT.XI)  SIGN=1.0 

IF(X.EQ.XI)  SIGN=0.0 

DO  50  N  =  l,  15 

DELMU=2*A/(6*N+3) 

AMU=0.0 

SUiM=0.0 

F0=(Q1( Y, ETA, XX, AMU) -Ql 0( Y , ETA , XX) )/( AMU-A) 

LL=6*N+3 

DO  40  NN=1,LL 

IF(ABS( AMU-A). LT.. 00001)  GO  TO  10 

IF(A3S( AMU+DELMU-A) .LT.. 00001)  GO  TO  10 

IF<A8S(  Av;U  +  DELMU/3.-A).LT..00001  )  GO  TO  10 

IF( ABSl AMU+2.*DELMU/3.-A) .LT..00001)  GO  TO  10 

Fl=(gi( Y, ETA, XX, AMU+DELMU/3. )-QlC(Y,ETA,XX) I /(AMU* 
1DELMU/3.-A) 

F2=(Q1 ( Y, ETA, XX, AMU+2 .*DELMU/3 . } -Q10 ( Y , ETA, XX ) ) / 
l(AMU+2.*DELMU/3.-A) 

F3=(Q1( Y, ETA, XX, AMU+OE LMU ) -Q10( Y , ETA ,XX ) )/( AMU+DELMU 
1-A) 

GO  TO  30 
10  F0=(Q1( Y, ETA, XX, AMU+DELMU } - Q10 ( Y , ETA , XX ) ) /{ AMU+DELMU 
1-A) 

GO  TO  40 
30  SUM=    (DELMU/3.  )*(F0 +3  .*FH-3.*F2+F3) +SUM 

F0=F3 
40  AMU=AMU+DELMU 

TEST(N) =SUM 

IF(N-l)  50,50,45 
45  MN=N-1 

MM=6*MN+3 

IF(ABS(  (TEST(iM)-TEST(N-l)  ) /TEST(N)  )  .LT.  .010)  GO  TO  60 
50  CONTINUE 
60  CONTINUE 

PV1F=2.*SUM 

EVALUATION  OF  THE  FINITE  INTEGRAL  IN  GN  USING 
2-A/MM  SUBDIVISION  SIZE 

DELMU=2*A/MM 

AMU=0.0 

SUMX=0.0 

SUMY=0.0 

F0=(Q2X(Yf ETA, X, XI ,AMU)-Q2  0X(Y ,ETA,X,XI ) )/{ AMU-A) 

Y0=(Q2Y(Y,ETA,X,XI  ,  AMU ) -Q2  0Y ( Y , ETA ,X , X  I  )  )/( AMU-A) 

DO  80  NN=1 ,MM 

IF(ABS< AMU-A). LT.. 00001)  GO  TO  70 

IF<ABS(  AMU+DELMU-A)  .LT.. 00001)  GO  TO  70 

IF(ABS( AMU+DELMU/3.-A) .LT. .00001 )  GO  TO  70 

IF(ABS( AMU+2.*0ELMU/3.-A) .LT..00001)  GJ  TO  70 

F1=(Q2X(Y, ETA,X,XI,AMU+DELMU/3. )-Q2QX(Y,ETA,X,XI))/ 
HAMU+DELMU/3.-A) 

F2=(Q2X(Y,ETA,X,XI , AMU+2. *DLLMU/3 . ) -Q20X ( Y , ETA , X , XI ) )/ 
1( AMU+2.*DELMU/3.-A) 

F3=(Q2X(Y,ETA,X,XI  ,  AMU+-DELMU  ) -Q2  OX  ( Y  ,  ET  A  ,  X  ,  XI  )  )/ 
1< AMU+DELMU-A ) 

Y1=(Q2Y(Y,ETA,X,XI , AMU+DE LMU/3 . )-Q20Y(Y ,ETA,X,XI))/ 
1  I AMU+DELMU/b.-A) 

Y2=(Q2Y(Y,ETA,X,XI , AMU+2 .*DELMU/3 . J-Q20Y ( Y , ET A,  X  ,  X  I )  ) / 
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l(AMU+2.*DELMU/3.-A) 

Y3=(Q2Y(Y,ETA,X,XI, AMU+DELMU)^Q20Y (Y,ETA ,X, XI ))/ 
1 ( AMU+DELMU-A) 

GO  TO  75 
70  CONTINUE 

F0=(Q2X(Y,ETA,X,XI,AMU+DELMU)-Q2CX(Y,ETA,X,XI ) )/ 
1( AMU+DELMU-A) 

Y0=(Q2V (Y, ETA,X,XI,AMU+DELMU)-Q20Y(Y,ETA,X,XI ) )/ 
1 ( AMU+DELMU-A) 

GO  TO  30 
75  CONTINUE 

SUMX=SUMX+    (DELMU/8. )*(F0+3.*F1+3.*F2+F3) 

SUMY=SUMY+  (DELMU/8.  )*  I  YO  +3  .  *Y  1+3.*Y2  +  Y3  ) 

F0  =  F3 

Y0  =  Y3 
80  AMU=AMU+DELMU 

PV2FX=2.*SUMX 

PV2FY=2.*SUMY 

EVALUATION  OF  THE  FINITE  INTEGRAL  IN  GYY  USING 
2*A/MM  SUBDIVISION  SIZE 

DELMU=2*A/MM 

AMU=0.0 

SUMYY=0.0 

YY0=(Q2YY(Y,ETA,X,XI  , AMU) -Q20Y Y ( Y ,ETA , X , X  I ) )/( AMU-A) 

DO  350  NiW  =  l,MM 

IF(ABS( AMU-A). LT. 0.00001)  GO  TO  320 

IF  (A3S( AMU  +  DELMU-A)  .LT. 0.00001  )  GO  TO  3  20 

IF  ( A3S( AMU+QELMU/3. -A) .LT. 0.00001)  GO  TO  320 

IF  (ABSCAMU+2. *DELMU/3. -A) .LT.  0.00001)  GO  TO  320 

YYl=(Q2YY(Y,ETA,X,XI,AMU+DELMU/3.)-Q20YY(Y,ETA,X,XI ) )/ 
1 (AMU+DELMU/3. -A) 

YY2=  (  Q2  YY  <  Y  ,  ET  A  ,  X  ,  XI ,  A MU  »-2  .  *DE  L  MU/3. )  -3  2CYY  I Y  ,  ETA  t  X , 
IX  I )  )/( AMU+2.*DELMU/3.-A) 

YY3=(  Q2YY(Y,ETA,X,XI , AMU  +  DELMU ) -Q20YY ( Y , ETA,X,XI ) )/ 
1 (AMU+DELMU-A) 

GO  TO  325 
320  CONTINUE 

YY0=(Q?YY< Y, ETA,X,XI,AMU+DELMU)-Q20YY(Y,ETA,X,XI) )/ 
1( AMU+DELMU-A) 

GO  TO  3  50 
3  25  CONTINUE 

SUMYY=SUMYY+( DELMU/8. ) *( YY0+3. *YYl+3 .-YY 2+Y Y3 ) 

YY0=YY3 
350  AMU=AMU+OELMU 

PV2FYY=2.*SUMYY 

EVALUATION  OF  THE  INFINITE  INTEGRAL  IN  G,  GN,  AND  GYY 
SIMULTANEOUSLY 

AMU=2*A 
DELMUO=DELMU 

F0=Q1( Y,ETA,XX,AMU)/( AMU-A) 

F0X=Q2X(Y,ETAfXTXI , AMU) /( AMU-A) 

F0Y=Q2Y(Y,ETA,X,XI ,AMU)/( AMU-A) 

F0YY=Q2YY( Y,ETA,X,XI ,AMU) /( AMU-A) 

SUM=0.0 

SUMX=0.0 

SUMY=0.0 

SUMYY=0.0 

DO  100  NN=1,200 

DO  95  N  =  l, 20 

F1=Q1(Y, ETA, XX, AMU+DELMU/3. ) /( AMU+DELMJ/ 3 .- A) 

F1X=Q2X(Y,ETA,X,XI , AMU+DELMU/3 . )  /  (  AMU  +  DELMU/.J  .-A) 

F1Y=Q2Y(Y,ETA,X, XI , AMU+DELMU/3. ) / ( AMU+DE LMU/3. -A) 

F  1  Y  Y  =  Q2  YY  (  Y  ,  L'  T  A  ,  X  ,  X  1  ,  AMU  +  J  E  L  MU  /  3  .  )  /  (  AMU  +  DE  L  MU  /  3  .  -  A  ) 

F2=Q1(Y, ETA,XX, AMU+2.*DELMU/3. ) / ( AMU+2 . * DELMU/ 3 . - A ) 


F2X=Q2X 

i-A) 
F2Y=Q2Y 

1-A) 
F2YY=Q2 

1-A) 
F3=Q1(Y 
F3X=Q2X 
F3Y=U2Y 
F3YY=Q2 
SUM  = 
SUMX= 
SUMY= 
SUMYY=S 
F0  =  F3 
F0X=F3X 
F0Y=F3Y 
F0YY=F3 
ASrt=EXP 
IF (ASM. 

95  AMU=AMU 
86    TESKNN 

TESTT(N 
TESTTT( 
TESTYY( 
IF(NN-1 

97  CONTINU 
IFtASM* 
IF(ABS( 
IF(ABS( 

1GO  TO  1 

98  IF(ABS( 
IF(ABS< 

IGO  TO  1 

99  IF(ABS( 
I F(ABS( 

IGO  TO  1 

93  IF(ABS( 

IF  (ABS 

IGO  TO  1 

96  CONTINU 
GO  TO  L 

100  CONTINU 

102  WRITE(6 

103  FORMAT ( 
105  CONTINU 

PV1I=2. 

PV2IX=2 

PV2IY=2 

PV2IYY= 

DELMU=. 

IF     (J.E 

AMU=0.0 

SUM=0.0 

SUMX=0. 

SUMY=0. 

SUMYY=0 

F 0=0.0 

FOX=FUN 

FOY=FUN 

FOYY=FU 

DO    200 

IF (AMU* 

DO    195 

F1X=FUN 

F2X=FUN 

F3X  =  FU:M 

F1Y=FUN 

F2Y=FUN 

F3Y=FUN 

F1YY=FU 

F2YY=FU 


(Y,ETA,X,XI,AMU+2.*DELMU/3. ) /( AMU  + 2. *DELMU/3 . 
(Y, ETA,X 
YY( Yt  ETA 


, XI ,AMU+2.*DELMU/3.)/(AMU+2.*DELMU/3. 
, X,XI, AMU*2.*DELMU/3.)/( AMUf 2.*DELMU/3, 


f ETA, XX , AMU*  DELMU)/( AMU+DELMJ-A) 

(Y,ETA,X,XI,AMU+DELMU)/( AMU+DEL/U-A) 
( Y, ETA,X,XI , AMU+DELMU)/( AMU+DELMU-A) 
YY( Y,ETA,X,XI, AMU+DELMU) / ( AMU+DE LMU- A) 
<DELMU/3.)*(F0+3.*FH-3.*F2+F3)+SJM 
<DELMU/8.)*(F0X  +  3.*FlX  +  3.*F2XfF3X)«-SUMX 
(DcLMU/3.)*(F0Y  +  3.*FiY  +  3.*F2Y«-F3Y)*SUMY 
UMYY+OELMU/8.  )*(F0YY+3.*F1YY+3.*F2YY+F3YY) 


YY 
(AMU*(ET 

LT.  .0001 

+DELMU 

)=SUM 

Nl-SUMX 

NN)=SUMY 

NNJ=SUMY 

)  100,10 

E 

LT. .0001 

SUM) .LT. 

(TEST(NN 

00 

SUMX) .LT 

(TESTT(M 

00 

SUMY) .LT 

(  i E  ST TT ( 

00 

SUMYYJ .L 

( ( TESTYY 

00 

p 

05 

E 

,103) 

3X34HINFINITE  INTEGRAL  DID  NOT  CONVERGE) 

E 

*SUM 

,*SUMX 

,*SUMY 

2.* SUMY Y 

Z/H 

Q.25)  GJ  TO  240 


A+Y) )/AMU 
)  GO  TO  86 


Y 
0,97 

0)  GO  TO  105 

.0001 >  GO  TO  98 

)-TEST(NN-i )) /TEST (NN)) .GT..  001) 

.  .0001)  GO  TO  99 

M)-TESTT(NN-i) )/TESTT(NMJ ) .GT..001  ) 

. .0001)  GO  TO  93 

T. 0.0001)  GO  TO  96 

(NN) -TESTY Yl NN- 1 ) ) /TESTY Y ( NN  ))  .GT . .001 ) 


.0 

3X(  Y, 
3Y(  Y, 
N3YY( 
NN  =  1, 
H.GT. 
N=l,2 
3X(  Y, 
3X(  Y, 
3X(  Y, 
3Y(  Y, 
3Y(  Y, 
3Y(  Y, 
N3YY( 
N3YY( 


ETA 

ETA 

Y,E 

ICO 

5.  ) 

0 

ETA 

ETA 

ETA 

ETA 

ETA 

ETA 

Y  ,  t: 

Y,E 


,X,XI,AMU) 
, X,XI ,AMU) 
TA,X,XI ,AMU) 

DELMU=.l 


tX, 
,  X  , 
tX, 
iX, 

f  A  , 

rX, 

TA, 


, AMU+DELMU/3. ) 

, AMU*2.*DELMU/3. ) 

, AMU+DELMU ) 

, AMU  +  DELMU/3.  ) 
, AMU+2.*DELMb/3. ) 
, AMU+DELMU) 


TA,X,XI  ,AM'J  +  2.*DELMU/3.  ) 


86 


195 
196 


199 


3YY=FU 
F(AMU+ 

1-FUN1 

2=FUN1 

3=FUN1 

0    TO    1 

ONTINU 

1=EXP( 

AMU+DE 

2=EXP( 

(AMU+2 

3=EXP( 

Y*ETA/ 

ONTINU 

UM= 

UMX  = 

UMY  = 

UMYY=S 

0=F3 

0X=F3X 

OY=F3Y 

0YY=F3 

SM=EXP 

F(ASM. 

MU=AMU 

ESTINN 

ESTT(N 

ESTTT( 

ESTYY( 

F(NN-1 

ONTINU 

F(ASM. 

F(ABS( 

F(ABS( 

0    TO    2 


F 

I 

F 

F 

F 

G 
120    C 

F 
1( 

F 
1* 

F 

1* 

130    C 

S 

S 

S 

S 

F 

F 

F 

F 

A 

I 

A 

T 

T 

T 

T 

I 

C 

I 

I 

I 
/lGr 
tOo     I  r  \  A  vj  ^  ( 

IF(A6S( 
1G0    TO    2 

207  IF(ABS( 
IF(ABS( 

1G0    TO    2 

204  IF(ASS( 
IF    UBS 

1GT. .001 

208  CONTINU 
GO    TO    2 

200    CONTINU 

WRITE(6 

202    FORMAT ( 

205  CONTINU 
GINF=-2 
GXINF=2 
GYINF=2 
GYYINF= 
GNSING= 
IF  (I  .E 
I F< l.EQ 
AIJJ-I- 
THETA1= 
AINJJ=I 
THETA2= 
AJNII=I 
THETA3= 
THETA=T 
IF(THET 
1F(THET 
IF(THET 
GSING=D 

1THETA/2 
2) *ALOS( 
3+(THETA 
4( 180. *5 


N3YY( 
DELMU 
(Y,ET 
(Y,ET 
( Y  ,  E  T 
30 
E 

-{  AMU 
LMU/3 
-(  AMU 
.-DEL 
-(  AMU 
COSH( 
E 

(OELM 

(DEL 

(DEL 

UMYY  + 


YfETAfXiXI t AMU+DELMU) 

.LT..001)     GO    TO     120 
A, XX, AMU+DELMU/3. ) 
A,XX, AMU+2. *DE LMU/3. ) 
A, XX, AMU+DELMU) 


+DELMU/3 .  }*H)*COS(  ( AMU  +  DELMU /3 .  )*XX)* 
. ) *Y*ETA/COSH( ( AMU  +  DELMU/3.  )*H) 
+2.*DELMU/3. )*H)*CGS(  ( AMU+2 . *DE LMU/3 .  )*XX) 
MU/3. )*Y*ETA/COSH( I AMU  +  2 . *DE LMU/3 . ) *H ) 
+DELMU)*H)*COS( (AMU+DE LMU)*XX)* (AMU+DELMU} 
( AMU+DELMU) *H) 

U/8. )*(F0+3.*F1+3.*F2+F3)+SUM 
MU/8.}*( F0X+3.*F1X+3.*F2X+F3X)+SUMX 
MU/8.)=MF0Y  +  3.*F1Y  +  3.*F2Y  +  F3Y)  +  SUMY 
(DELMU/8.)*(F0YY+3.*FIYY+3.*F2YY+F3YY) 


YY 

( AMU*(ETA+Y) ) 

LT..0001)     GO    TO    196 

+DELMU 

)=SUM 

N)=SUMX 

NNJ=SUMY 

NNJ=SUMYY 

.     200,200,199 

E 

LT. .00010)     GO    TO    205 

SUM) .LT. .0001  )     GO    TO    206 

( TEST (NM)- TEST (NN-1) ) /TEST ( NN ) ) . GT. . 00 10 ) 

00 

C  I  I  •  A  V  1        IT  r%  ,-^  n,  l    y        /~n      T  r~>      O  ^  "7 

Jdl/WaLteavJWw/lj        O  v_>        <    <_>      C  \J   I 

( TESTT(NN)-T£$TT(NN-1) ) / TESTT ( NN ) ) . GT. .00 10 ) 

00 

SUMY) .LT. .0001 )     GO    TO    204 

(  TESTTT(NN)-TESTTT(NN-1)  } /TESTTT  (  NiM  )  ).GT..0010) 

00 

SUMYY) .LT. 0.0001)     GO    TO    203 

(  (TESTY Y(  'JN)-TESTYY(NN-l)  )/(TESTYY(NN))). 

)    GO    TO    200 

E 

05 

E 

,202) 

3X14HNG    CONVERGENCE) 

E 

*SUM 

•*SUMX 

,*SUMY 

2.*SUMYY 


Q. 

.J 
J 

A3 
+  N 
A3 
-J 
A3 
HE 
A2 
A3 
A. 
EL 
.  ) 
TH 
/2 
.) 


25)    GO    TO    218 
}    GO    TO    220 


S(  AIJ 

PTS-J 
S  (  A  I N 
-NPTS 
S(  AJN 
TA1 
.LT 
.LT 
GT. 
THE* 
*ALOG 
ETA/ 2 
. -DEL 
+  (THE 


J)*DELTHE 

JJ  )*DELTHE 

I  I  )*DELTHE 

HETA)     THETA=THETA2 

HETA)     THETA=THETA3 

5)     GO    TO    218 

LUG(  2.  ) +  2.* (-DEL  THE/ 2.  + (DEL  THE/ 4.+ 

( THETA/2.+DELTHE/4.)-( THETA/2.-DELTHE/4. 

. - D ELT  H  E  /  4  .  )  -  (  0  E  L  T  H  E  /  <*  .  +  T  H ETA/2. )  *  *  3/18. 

THE/4. ) **3/18.~< THETA/2.+DELTHE/4. )**5/ 

TA/2.-DELTHE/4. )**5/( 180.* 5. )-(THETA/2.+ 


5DEL THE/4. ) 

6(2635. *7.) 
GXSING=(X- 
GYSING=(Y- 
GO  TO  230 
218  GSING=DELT 
GXSING=(X- 
GYSING= (Y- 
GO  TO  2^>0 
220    CONTINUE 

GSING=OELT 

1( 18.*16.)- 

2(2835. *7.* 
GXSING=GNS 
GYSING=GNS 
230  CONTINUE 

G   =GSI NG/ 

H-PV1I+GINF 
GX=GXSING- 

1FUN4(Y ,ETA 
GY=GYSING- 

IFUN4(Y,ETA 
IF  (  I  .  L  E  .  N 
FURYY=1 ,/S 

1( Y-ETA) )/S 
FUNYY=1./S 

MY  +  ETA)  )/S 
GYY=FURYY+ 

1FUN2( Y ,ETA 
235  CONTINUE 

IF  (I.EQ.2 
GN=GNSING- 

l*ANX+( PV2F 

2SIGN)#CMPL 
GO  TO  250 
240  CONTINUE 

G=-(PV1F+P 
GN  =  -( PV2FX 

lANXf ANY tSI 
2  50  CCNT1 NUE 
RETURN 
END 


**7/( 2  835.*7.)+( THETA/2.-DELTHE/4.)**7/ 

) 

XI)/((X-XI)**2+(Y-ETA)**2) 

ETA) /(( X-XI )**2+( Y-ETA )**2) 

HE*ALOG(SQRT( ( X-X I ) **2  +  ( Y-ETA ) **2  )  ) 
XI )/(<X-XI)**2+(Y-ETA)**2) 
ETA)/( (X-XI )**2+( Y-ETA )**2) 


HE*AL0G(DELTriE/2. ) -D EL THE- ( DE L 

(DELTHE**5)/(1S0.*5.*2  56.)-(DE 

2  56.- lb.  ) 

INb*ANX 

ING*ANY 


THE**3)/ 
LTHE**7)/ 


DcLTH 
+CMPL 
FUNP( 
tXX,l 
FUNR( 
,XX,0 
PTS) 
WRT(  ( 
QRT(  ( 
ORT(  ( 
QRT(  ( 
rUNYY 
,  XX  )  -• 


E-A 
X(0 
X,Y 
.0, 
X,Y 
.Or 


LOG( 
.,-1 

tXI  , 
0.0, 
»XI  , 
3, 


1. 

TO  2 
XI  )* 
XI  )* 


(X- 
(X- 
(X-XI )* 

(X- 
+  PV 
A*A 


XI  ) 
'2FYY 


SQRT( 
.  )*FU 
ETA,1 
SIGN) 
ETA,0 
SIGN) 
35 

*2  +  (Y 
*  2  *-  (  Y 
*2f  (  Y 
*2+(Y 
+  PV2I 


(X-XI)**2+(Y+ETA)**2) )+PVlF 

N2( Y, ETA, XX) 
.0,0.0)+PV2FX 
*CMPLX(0.0,-1 
.Otl .0)+PV2FY 
*CMPLX (0.0,-1 


fPV2IXfGXINF+ 

.0) 

+PV2IY+GYINF* 

.0) 


-ETA)**2)**3) 
-ETA)**2)**5) 
+ETA)**2)**3 ) 
+ETA)**2)#*5) 

YY+GYYINF+CW 


5)  GO  TO  250 

FJNR{X,Y,XI » ETAfANX,ANY)+(PV2F 
Y+PV2IY  +  GYINF)*ANY4-FUN4(Y,ETA, 
X  (  0  .  0  ,  -  1 .  0  ) 


VII )+CMPLX(0. t-1. )*FUN2(Y,ETA, 

+  PV2IX)*ANX-(PV2FY  +  PV2IY)*A,NY  + 
GN)*CMPLX<0. 0,-1.0) 


-(( Y-ETA)**2+ 
-( ( Y+ETA)**2+ 
LX(0. 0,-1.0)* 


X+PV2IX+GXIMF) 
XX, ANX, ANY, 


XX) 

FUN4( Y, ETA, XX, 


THIS  SUBROUTINE  CALCULATES  G  AND  ITS  DERIVATIVES, 
GXfGYfGNi  AND  3YY,  BY  USE  OF  THE  SERIES  FORM  FOR  THE 
CASE  (X(I)  -  X(J))  GREATER  THAN  SMIN 

SUBROUTINE  GREENS  ( A , ANU, SH2 AH, SHAH ,CHAH ,CGSAMU, S INAMU 
1 , AMU,COEFG,SHYI , CHYI, SHYJ,CHYJ, I , J,XV,VV,XC,YC, AN  XI , 
2ANYIf  ANXJi  ANYJiGiJ,GNI  J,G.\IJI,aXIJ»GXJI,GYI  J,GYJI  ,GYY) 

COMPLEX  GUf  GNIJ,GNJI  ,GXI  J,GXJ  I  ,GYIJ,GYJI  ,  GYY 

DIMENSION  COS A MU( 200,2  5) , SINAMUt  200,25)  , AMU (200)  , 
1C0EFG(20C) 

DIMENSION!  TESTK40)  ,TEST2(40)  ,TEST3(40),TEST30(40), 
1TEST4(40) 

COMMON/CONST/H, J , DE LTHE , SMI N ,NPT S , NSPTS 

CCMMON/CPI/PI 

SUM1=0.0 

SUM2=0.0 

SUM3=0. 0 

SUM30=0.0 

SUM4=0.0 

DO  20  N=l,40 

DO  15  KK=1,5 

K=(N-1) *5  +  KK 

SNI=SINAMU(K,I ) 

SNJ=SINAMU( K, J ) 

CSI=COSAMU(K,I) 

CbJ=COSAMU(K, J) 
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AK=AMU( K) 

CCK=COEFG(K) 

VAL=AK*ABS(XV-XC) 

IF    (VAL. GE. 75.0)     GO    TO    10 

EXIJ=EXP(-VAL) 

GO    TO    1 1 

10  EXIJ=0.0 

11  CONTINUE 

SUM1  =  SUMH-C0K*CSI*CSJ*EXIJ 
SUM2=SUM2+CQK*CS1 *CSJ*EXI J* AK 
SUM3=SUM3+COK*AK*SNI*CSJ*EX  IJ 
SUM30=SUM30+C0K*AK*SNJ*CSI*EXIJ 
SUM4=SUM4+AK*AK*C0K*CSI*CSJ#EXI J 

15  CONTINUE 
TEST1(N)=SUM1 
TEST2(  N)  =  SU,M2 
TEST3(N)=SUM3 
TEST30(N)=SUM30 
TEST4(N)=SUM4 
IF(N-l)     20,20,16 

16  IF(ABS(  SUM1) .LE. 0.000001)     30    TO    17 
IFUBSUTESTl(i\l)-TESTl(N-l>)/TESTKN)).GT..0010)     GO 

1T0    20 

17  IF(ABS( SUM30J.LE .0.000001)     GO    TO    18 

IF(ABS( (TEST30(N)-TEST30(N-1) ) /TESTS 01 N) ) .GT. .0010) 
1G0    TO    20 

18  IF(ABS( SUM3) .LE. 0.000001)     GO    TO     19 

IF(ABS( (TEST3(N)-TEST31N-1) )/TEST3(N) J.GT..0010)     GO 
1T0    20 

19  IF(ABS(  SUM2) .LE. 0.000001)     GO    TO    21 

IF(ABS(  (  TEST2(N)-TEST2(.\-1)  )/TEST2(N))  .GT..0010)     GO 
1T0    20 
21     IF(ABS( SUM4) .LE. 0.000001)     GO    TC    30 

IF(ABS( <  TEST4<N)-TEST4(N-1) )/TEST4(N) ) .GT..0010)     GO 
1T0    20 

GO    TO    30 

20  CONTINUE 
WRITE(6,25)     I, J 

25    FORMATQ0X23HGREENS    DID    NOT    CONVERGE  ,  2X2HI  =  I  2  , 

12X2HJ=I2) 
30    CONTINUE 

IF     (XV.LT.XC)     SIGN=-l.O 

IF     (XV.GT.XC)     SIGN=1.0 

IF    (XV.EQ.XC)     SIGN=0.0 

TERM4=4.*PI*CHYI*CHYJ*SIN( A*A6S(XV-XC) ) / ( 2. *A*H* SH2AH) 

TERM5=4.*PI*CHYI*CHYJ*CCS( A*ABS< XV-XC) )/ ( 2. *A*H+SH2AH) 

TERM6=A*TERM5*SIGN 

TERM7=A*TERM4*SIGN 

TERM8=4.*PI*A*SrtYI*CHYJ*C0S(A*ABS(XV-XC ) ) /<2.*A*H+ 
1SH2AH) 

TERM9=4.*PI*A*SHYI*CHYJ*SIN(A*ABS(XV-XC) )/<  2.*A*H  + 
1SH2AH) 

TERM8Q=4.*PI*A*SHYJ*CHYI*C0S(A*ABS(XV-XC ) J/(2.*A*Hf 
1SH2AH) 

TERM90  =  4.*PI*A*SHYJ=::CHYI*SIN(A*A&S{XV-XC  ) )/(2.*A*H+ 
1SH2AH) 

TERM10=A*A*TERM4 

TERM11=A*A*TERM5 

IF  (J.EQ.25)  GO  TO  70 

G IJ=C  MD  LX (  SUM1  +  T  E  RM4 ,  -T  ER M5 ) 

GXIJ=    CMPLX(-SIGN*SUM2fTERM6fTERM7) 

GYIJ=CMPLX(-SUM3+TERM9,-TtRM8) 

IF  ( I.EQ.25)  GO  TO  40 

GX JI=C  MPLX ( S IGN*  SUM2-TE RM6 ,-TERM7 J 

GYJI=CMPLX(-SUM30+TERM90,-TERM80) 

GNI  J  =  CrtPLX(-ANXI*SIGN*SUM2  +  ANXI*TERr16-ANYI*SUM3*-ANYI* 
1TERM9,  ANXI*TERM7-ANYI*TERM8 ) 

GNJI=CMPLX(ANXJ*SIGN*SUM2-ANXJ*TERM5-ANYJ*SUM30+ANYJ* 
lTERM90t-ANXJ*TERM7-ANYJ*TERM80) 

GO  TO  60 
40  CONTINUE 

GYY=CMPLX(-SUM4+TERM10|-TERM11 ) 


89 


60  CONTINUE 

GO  TO  SO 
70  GIJ=CMPLX(-SUM1-TERM4,-TERM5) 

GNIJ=CMPLX(ANXi*SIGN*SUM2-ANXI*TERM6+ANYI*SUM3-ANYI* 
LTERM9, ANXI*TERM7-ANYI*TERM8) 
80  CONTINUE 

RETURN 

END 

THIS  SUBROUTINE  INVERTS  COMPLEX  MATRICES  TO  SOLVE 
THE  MATRIX  EQUATION  ALPHAi I , J ) *F ( I , 1 )  =  HH(I,1)  AND 
ALPHAt I  ,J)*F1( 1,1)  =  PK(I ,1) 

SUBROUTINE  COMAT ( N ,M , A , 6, Q , I ) 

INTEGER  C,H,R,Q,Z 

COMPLEX  A,3,0,TT,P 

DIMENSION  A(N»N) ,6(N,M) ,C( 100,3) 

D  =  ( 1  .0,0.0) 

DO  20  J  =  i,N 

20  C(J,3)  =  0 

DO  21  K  =  1,N 

TT  =  (0.0,0.0) 

T  =  0.0 

DC  4  J  =  1,N 

IF  (C( J, 3)  .EQ.  1)  GO  TO  4 

DO  5  H  =  i,N 

IF  (C(H,3)  ~  1)  15,5, 12 

15  IF  (T  .GE.       CABS( A( J,H) ))   GO  TO  5 
R  =  J 

Q  =  H 

T  =       CABS( A( J,H) ) 
5  CONTINUE 
4  CONTINUE 

C(Q,3)  =  C<Q*3)  +  1 

C(K,1)  =  R 

C(K,2)  =  Q 

IF  (R  . EQ.  Q)  GO  TO  11 

D  =  -D 

DO  6  L  =  1  ,N 

TT  =  A(R,L) 

A(R,L)  =  A(Q,L) 
8  A(Q,L)  =  TT 

IF(M  .LE.  0)  GO  TO  11 

DO  2  L  =  1 ,M 

TT  =  B( R,L) 

B(R,L)  =  B(Q,L) 

2  B(QtL)  =  TT 
11  P  =  A(Q,Q) 

A(Q,Q)  =  (1.0,0.0) 
DO  13  L  =  1,N 
13  A(Q,L)  =  A(0,L)/P 

IF  (M  .LE.  0)  30  TO  1 
DO  3  L  =  1,M 

3  B(Q,L)  =  B(Q,L)/P 
1  DO  21  Z  =  1,N 

IF  (Z  . EQ.  Q)  GO  TO  21 
TT  =  A( Z,Q) 
A(Z,Q)=(0. 0,0.0) 
DO  16  L  =  1,N 

16  A(Z,L)  =  A(ZTL)  -  A(Q,L)*TT 
IF  (M  . LE.  0)  GO  TO  21 

DO  17  L  =  1,M 

17  B(ZtL)  =  B(Z,L)  -  B(Q,L)*TT 

21  CONTINUE 

DO  19  II  =  1,N 

L  =  N  +  1  -  II 

IF(C(L,1)  .EQ.  C(L,2))  GO  TO  19 

R  =  C(L,1) 

Q  =  C(L,2) 
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DO  7  K  = 

ItN 

TT  =  A<  K, 

rR) 

A(KfR)  = 

A(K,  Q 

7 

AlKtQ)  = 

TT 

19 

CONTINUE 

DO  18  K  = 

=  L.N 

IF  (C(K,3)  .NE 

18 

CONTINUE 
I  =  1 

50 

RETURN 

12 

I  =  2 

GO  TO  50 
END 

//GO.: 

SYSIN    DD  * 

1)    GO    TO    12 
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